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Preface 


The  First  North  American  Workshop  on  Modeling  the  Mechanics  of 
Off-Road  Mobility  was  held  5  and  6  May  1994  at  the  U.S.  Army  Engineer 
Waterways  Experiment  Station  (WES)  in  Vicksburg,  Mississippi.  The 
workshop  was  sponsored  by  the  U.S.  Army  Research  Office  under  the 
Terrestrial  Science  Program  of  the  Engineering  and  Environmental  Sciences 
Division. 

The  idea  for  this  workshop  originated  with  Mr.  David  A.  Homer,  Mobility 
Systems  Division  (MSD),  Geotechnical  Laboratory  (GL).  The  workshop  was 
subsequently  organized  by  Mr.  Roger  W.  Meier,  MSD,  under  the  general 
supervision  of  Mr.  Newell  R.  Murphy,  Chief,  MSD,  and  Dr.  William  F. 
Marcuson,  III,  Chief,  GL.  This  report,  which  documents  the  proceedings  of 
the  workshop,  was  edited  by  Messrs.  Meier  and  Homer  with  the  assistance  of 
Mr.  Jody  Priddy,  MSD.  Mr.  Meier  wrote  the  introduction  and  workshop 
summary  in  Chapter  1 .  Ms.  Sally  Shoop,  U.S.  Army  Cold  Regions  Research 
and  Engineering  Laboratory,  provided  the  synopsis  of  the  plenary  discussions 
in  Chapter  2. 

The  workshop  organizers  wish  to  thank  Drs.  Peter  Haff,  Duke  University, 
Shrini  Upadhyaya,  University  of  California  at  Davis,  and  Paul  Corcoran, 
Caterpillar,  Inc.,  for  serving  as  working  group  chairmen  and  providing  us  with 
synopses  of  their  group’s  discussions. 

Dr.  Robert  W.  Whalin  was  Director  of  WES  during  the  preparation  and 
publication  of  this  report.  COL  Bruce  K.  Howard,  EN,  was  Commander. 


The  content!  of  this  report  are  not  to  be  used  for  advertising,  publication, 
or  promotional  purposes.  Cite,. ion  of  trade  names  does  not  constitute  an 
official  endorsement  or  approval  of  the  use  of  such  commercial  products. 


1  Introduction 


Background 

In  the  current  atmosphere  of  belt-tightening  and  streamlining,  the  Army 
needs  a  means  of  assessing  the  mobility  performance  capabilities  of  candidate 
next-generation  vehicles  without  first  spending  tens  of  millions  of  dollars 
building  vehicle  prototypes.  Shortcomings  and  design  flaws  that  compromise 
mobility  performance  must  be  identified  before  they  are  incorporated  into  the 
prototypes  at  Army  expense. 

The  task  of  evaluating  the  mobility  performance  of  vehicles  in  currently 
accomplished  in  part  with  the  NATO  Reference  Mobility  Model  (NRMM).1 
NRMM  uses  a  collection  of  algorithms,  numerics,  and  empirical  relationships 
to  forecast  maximum  steady-state  vehicle  speed  as  a  function  of  driver, 
vehicle,  weather,  and  terrain  characteristics.  The  empirical  relationships 
embody  more  than  40  years  of  vehicle  mobility  research,  testing,  and 
evaluation.  Efforts  to  improve  NRMM  and  expand  its  empirical  database 
continue  to  this  day. 

There  are,  however,  advances  in  vehicle  design  on  the  horizon  that  will 
produce  vehicles  that  operate  outside  the  limits  of  NRMM’s  empirical 
database.  Many  of  these  concepts — such  as  central  tire  inflation  systems, 
rubber  belt  traction  elements,  active  suspensions,  and  weight-saving  reactive 
armor — are  already  appearing  in  production  and  demonstration  vehicles.  Other 
concepts — such  as  lightweight-composites  and  appliqud  armor,  zero-ground 
pressure  running  gear,  and  electric  drive  technologies — are  just  around  the 
comer.2  In  many  cases,  the  performance  characteristics  of  these  vehicles 
simply  cannot  be  described  within  the  bounds  of  the  existing  NRMM  database. 
Expanding  NRMM’s  empirical  database  to  include  these  advanced  vehicles 
would  involve  extensive  mobility  field  tests  on  prototype  vehicles.  This  is 
exactly  the  type  of  expense  that  the  Army  can  no  longer  afford. 


1  Richard  Ahlvin  and  Peter  Haley.  1992.  “NATO  Mobility  Model  Edition  II.  NRMM  II  User's  Guide.” 
Technical  Report  GL-92-19.  U.S.  Army  Engineer  Waterways  Experiment  Station.  Vicksburg.  MS. 

3  Board  on  Army  Science  and  Technology.  1992.  “STAR  21:  Strategic  Technologies  for  the  Army  of  the 
Twenty-First  Century.  Mobility  Systems.”  National  Academy  Press.  Washington.  DC. 
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To  accurately  access  future  vehicle  performance  characteristics  in  the 
absence  of  costly  prototyping  and  field  testing,  the  Army  must  develop  an 
ability  to  perform  virtual  mobility  testing — the  simulation  of  vehicle  mobility 
performance  in  the  computer.  This  will  require  advances  in  the  numerical 
modeling  of  vehicle  dynamics  and  vehicle-terrain  interaction  and  in  the 
characterization  of  the  terrain  for  modeling  purposes. 

In  order  to  1)  assess  the  current  state  of  the  art  in  vehicle  mobility 
modeling,  2)  identify  the  most  promising  areas  of  current  research,  and  3) 
determine  the  most  profitable  directions  for  future  research,  the  Mobility 
Systems  Division  of  the  U.S.  Army  Engineer  Waterways  Experiment  Station 
(WES)  invited  recognized  leaders  in  the  field  of  vehicle  mobility  modeling 
from  throughout  the  United  States  and  Canada  to  participate  in  a  two-day 
workshop  on  “Modeling  the  Mechanics  of  Off-Road  Mobility.”  This  report 
documents  the  proceedings  of  that  workshop. 


Workshop  Summary 

Prior  to  the  workshop,  participants  were  asked  to  provide  a  brief  (3-4  page) 
technical  note  or  extended  abstract  describing  their  current  mobility  modeling 
research  efforts.  These  submissions,  which  are  reproduced  in  this  report,  were 
assembled  into  a  workshop  preprint  volume  that  was  provided  to  all  of  the 
participants  when  they  arrived  at  the  workshop  site.  The  technical  notes  were 
used  by  the  workshop  organizer  to  determine  the  interests  of  each  participant 
in  order  to  assign  them  to  individual  working  groups.  The  preprint  volume 
also  served  to  “introduce”  participants  to  one  another.  Because  the  participants 
came  from  a  wide  variety  of  organizations  and  had  a  wide  variety  of 
backgrounds,  not  all  of  them  knew  each  other  or  were  familiar  with  each 
other’s  research  work.  It  was  hoped  that  the  preprint  would  help  initiate  off¬ 
line  discussions  between  researchers  that  might  lead  to  fruitful  collaborations. 

The  workshop  began  with  a  brief  welcome  by  the  Commander  and  the 
Director  of  WES  and  opening  remarks  by  the  workshop  organizer.  These  were 
followed  by  two  keynote  speeches,  which  are  included  here.  The  first  speech, 
by  Prof.  J.  Y.  Wong  from  Carleton  University  in  Ontario,  Canada,  addressed 
the  state-of-the-art  in  the  analytical  modeling  of  steady-state  off-road  mobility 
for  wheeled  and  tracked  vehicles.  The  second,  by  Dr.  Shrini  Upadhyaya  from 
the  University  of  California  at  Davis,  described  his  recent  research  into  the 
backcalculation  of  in  situ  soil  properties  from  field  test  results  using  response 
surface  methodologies. 

The  afternoon  was  spent  in  a  group  discussion  trying  to  determine  which  of 
the  existing  mobility  modeling  paradigms  did  and  did  not  work  so  we  could 
better  determine  the  areas  where  additional  research  was  needed.  A  synopsis 
of  that  session,  presented  in  the  next  chapter,  has  been  provided  by  the  session 
chairperson,  Ms.  Sally  Shoop  from  the  U.S.  Army  Cold  Regions  Research  and 
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Engineering  Laboratory  (CRREL).  During  that  group  discussion,  the  lack  of 
any  standardization  of  in  situ  tests  for  ascertaining  vehicle  mobility  was  noted. 
A  small  working  group  broke  off  from  the  main  group  to  address  that  issue. 
Mr.  George  Mason  from  WES  has  submitted  a  synopsis  of  those  discussions. 
That  synopsis  is  paraphrased  here: 

Any  standardized  test  for  determining  in  situ  soil  properties  must 
account  for  1 )  soil  layering  and  heterogeneity,  2)  the  effects  of  changes 
in  moisture  content  over  depth,  3)  the  loss  of  strength  that  results  from 
remolding,  4)  the  effects  of  organic  matter  such  as  grass  and  roots,  and 
5)  the  critical  depth  at  which  soil  shear  properties  will  dictate  mobility. 

It  must  also  permit  the  backcalculation  of  unique  soil  properties,  unlike 
tests  such  as  the  military  cone  penetration  test  which  can  produce  the 
same  test  results  for  many  different  combinations  of  soil  strengths  and 
compressibilities. 

A  standardized  test  should  1)  be  able  to  measure  changes  in  strength 
with  depth  at  5  cm  increments,  2)  produce  shearing  patterns  that 
correlate  in  some  way  with  the  width  of  loaded  area  of  interest  (e.g.,  of 
a  track  pad  or  a  tire  contact  patch),  3)  include  a  definition  of  the 
remoldability  of  the  soil,  4)  be  capable  of  rapid  measurements  in  low 
soil  strengths,  and  5)  be  transportable  to  the  field. 

The  plenary  session  was  concluded  by  recommending  that  the  remainder  of 
the  workshop  be  devoted  to  the  identification  of  specific  research  needs  and 
the  facilitation  of  possible  cooperative  research  between  the  participants. 

To  that  end,  the  participants  were  broken  up  into  three  working  groups  of 
equal  size  the  next  morning.  Each  group  was  asked  to  return  two  hours  later 
with  an  enumerated  list  of  the  five  biggest  knowledge  gaps  in  vehicle  mobility 
modeling.  They  were  also  asked  to  address  the  question:  “If  those  gaps  were 
to  be  filled  tomorrow,  what  would  it  buy  us?”  TTie  participants  were  assigned 
to  the  different  working  groups  based  on  perceived  similarities  in  their  needs 
and  interests.  For  example,  one  group  contained  most  of  the  researchers 
mainly  interested  in  ride  dynamics.  Another  group  was  composed  of  the 
researchers  primarily  concerned  with  the  mechanics  of  vehicle/terrain 
interaction. 

In  all,  17  different  knowledge  gaps  were  identified  (two  of  the  groups 
submitted  lists  with  more  than  five  knowledge  gaps  and  several  gaps  were 
identified  by  more  than  one  group): 

•  need  to  determine  where  the  existing  modeling  techniques  don't  work 

•  need  to  relate  soil  mechanics  properties  to  intrinsic  oil  properties 

•  need  data  on  the  spatial  and  temporal  variability  of  soils 
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•  need  better  ways  of  coping  with  inhomogeneity 

•  need  better  ways  to  measure  and  evaluate  soil  properties  in  situ 

•  i»<  .d  valid,  repeatable  (standard)  characterization  of  the  real  world 

•  need  to  instrument  the  soil  without  affecting  its  structure 

•  need  to  measure  and  understand  dynamic  soil  properties 

•  need  to  measure  soil  adhesion  and  understand  its  role  in  mobility 

•  need  to  estimate  changes  in  soil  resulting  from  vehicle  traffic 

•  need  an  adequate  effective  stress  theory  for  multi-phase  media 

•  need  to  understand  the  mechanics  of  continuous  (large  strain)  failure 

•  need  to  better  understand  the  correlation  between  slip  and  sinkage 

•  need  to  predict  interface  stresses  during  acceleration  and  braking 

•  need  to  determine  the  3-D  response  of  terrain  to  vehicle  steering 

•  need  a  good  model  for  RAMD  prediction 

•  need  numerical  models  or  lookup  tables  to  speed  model  execution 

If  those  knowledge  gaps  were  filled,  the  following  could  be  accomplished: 

•  real  time  simulation  of  vehicle  mobility 

•  validation  of  existing  and  new  mechanisms  and  models 

•  accurate  modeling  of  vehicle  agility/maneuverability 

•  the  ability  to  predict,  design  for,  and  control  soil  compaction 

•  the  ability  to  predict  the  effects  of  tread  and  grouser  design 

•  prediction  of  site-specific  vehicle  mobility  (e.g.,  virtual  test  course) 

•  prediction  of  vehicle  failure  and  estimates  of  reliability 

•  enhanced  mobility-based  design  and  procurement 

There  are  far  fewer  “knowledge  uses”  than  there  are  knowledge  gaps 
because  many  of  the  participants  had  very  similar  needs  and  desires  and  some 
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of  those  needs  could  only  be  met  if  several  knowledge  gaps  were  filled.  It  is 
somewhat  surprising  that  there  were  almost  as  many  knowledge  gaps  (17)  as 
there  were  participants  (22)  despite  the  similarity  in  their  knowledge  needs. 
Perhaps  that  is  a  strong  indication  that  there  is  still  much  work  to  be  done  in 
understanding  and  modeling  the  mechanics  of  off-road  mobility. 

Despite  the  lack  of  consensus  as  to  our  most  pressing  research  needs,  the 
workshop  was  successful  in  that  it  served  to  let  everyone  working  in  the  field 
know  where  they  stand  as  a  group.  The  current  state  of  the  art  was  illustrated, 
ongoing  research  was  brought  to  light  and  discussed,  and  knowledge  gap;  that 
need  to  be  filled  through  future  research  were  identified.  This  was  especially 
important  for  the  several  participants,  invited  at  the  behest  of  ARO,  whose 
backgrounds  were  outside  the  realm  of  vehicle  mobility  modeling.  Those 
researchers  proposed  some  inventive  new  approaches  to  mobility  modeling. 
Several  of  those  approaches  are  currently  being  investigated  and  will  probably 
lead  to  research  contracts  or  cooperative  research  agreements.  If  nothing  else 
were  to  come  of  this  workshop,  that  alone  will  have  made  it  worthwhile. 

At  the  end  of  the  workshop,  there  was  unanimous  agreement  that  the 
workshop  should  become  a  recurring  event.  ARO  has  agreed,  in  principle,  to 
sponsor  a  recurring  workshop.  The  Mobility  Systems  Division  of  WES  has 
agreed  to  host  it  again.  That  “Second  North  American  Woricshop  on  Modeling 
the  Mechanics  of  Off-Road  Mobility”  is  tentatively  scheduled  for  August  1995. 
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2  Plenary 


This  plenary  session  was  loosely  structured  with  the  broad  objectives  of 
concentrating  on  identifying  the  knowledge  gaps  in  the  current  state  of 
modeling  vehicle-terrain  interaction.  Because  we  were  of  very  diverse 
backgrounds,  and  with  a  wide  variety  of  applications  for  such  modeling,  the 
discussion  began  by  grouping  the  applications  as  those  dealing  with  prediction 
and  improvement  of  performance  (such  as  vehicle,  tire,  track,  compactor,  and 
agricultural  tools),  and  those  dealing  with  the  resulting  soil  deformation  (such 
as  compaction,  mass  flow,  and  terrain  damage). 

The  discussion  then  focused  on  the  general  needs  of  vehicle-terrain  models, 
as  follows: 

•  The  need  to  relate  remotely-sensed  mapped  parameters  to  soil 
parameters  and  vehicle  performance  models. 

•  The  need  to  use  FEM  (or  other  sophisticated  tools)  to  generate  lookup 
tables  for  bigger,  virtual  reality  type  models  which  must  run  in  real 
time. 

•  The  need  for  a  common  or  standard  measurement  device,  and/or  a 
“standard”  piece  of  ground  to  validate  soil  devices  and  models.  The 
general  consensus  was  for  coming  up  with  a  “standard”  shear  device. 

•  It  was  generally  agreed  that  although  the  vehicle  input  needed  for  such 
models  is  well  defined  and  understood,  the  soils  characterization  and 
modeling  has  a  long  way  to  go. 

At  this  point,  a  small  group  split  off  to  discuss  the  device  standardization 
issue.  The  remaining  participants  discussed  the  following  issues  regarding  the 
specific  capabilities  that  are  lacking  to  do  the  above.  Some  of  this  discussion 
continued  the  second  d?y. 

•  The  3-D  contact  between  the  driving  element  and  the  soil  is  not  often 
known  but  must  be  used  as  input  for  rigorous  numerical  modeling. 
(However,  using  contact  mechanics,  it  is  the  material  properties  that 
must  be  well  known  and  the  contact  is  calculated.) 
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•  Is  it  valid  to  apply  the  traditional  Janosi's  shear  equation  and  Bekker’s 
sinkage  equation  for  the  case  of  a  3-D  soil-tire  contact  surface? 

Janosi's  equation  assumes  horizontal  shear  and  Bekker's  sinkage 
assumes  vertical  (or  normal  to  the  loaded  surface)  deformation. 
However,  the  stresses  and  deformations  involved  in  soil-tire  contact  on 
deformable  terrain  are  not  necessarily  vertical  and  horizontal.  As 
Bekker's  equation  assumes  that  the  pressure  is  hydrostatic  the  direction 
is  not  critical,  however,  the  asymmetry  of  the  pressure  distribution 
beneath  the  wheel  causes  problems.  Janosi's  equation  is  based  on  the 
shear  and  normal  forces  on  the  wheel  being  vertical  and  tangential  and 
was  intended  to  be  used  with  the  sum  of  the  forces  on  the  wheel. 

•  Soil  energy  absorption  terms  needed  for  vehicle  dynamics  simulations 
are  lacking. 

•  A  good  description  and  understanding  of  slip-sinkage  is  needed. 

•  Description,  numerical  formulation  and  modeling  methodology  for 
dealing  with  the  rate  effects  associated  with  dynamic  loading  of  very 
wet  and  saturated  soils. 

•  Problems  associated  with  modeling  large  deformations  are  encountered 
in  FEM.  The  use  of  Eulerian  FEM  codes  was  discussed,  as  were 
Discrete  Elements  and  particle  theory  along  with  their  large 
computational  requirements. 

•  The  ability  to  describe  and  model  heterogeneous  material  (such  as 
boulders  and  roots)  is  needed. 

•  The  ability  to  handle  the  effects  of  multiple  vehicles,  loading  and 
unloading,  repetitive  loads  is  needed. 

•  A  good  (constitutive)  model  for  soil  under  dynamic  loads  is  needed. 

•  Money  is  needed. 

A  good  deal  of  the  remaining  discussion  centered  around  tight  budgets  and 
where  to  get  the  funding  to  develop  and  implement  these  projects.  Some 
funding  is  available  from  the  Army  Research  Office  (ARO)  and  some  is 
available  from  government  laboratories  through  their  respective  Broad  Agency 
Announcements  (BAA).  Other  sources  mentioned  were  ARPA  (Advanced 
Research  Projects  Agency),  AFOSR  (Air  Force  Office  of  Scientific  Research), 
and  CPAR  (Construction  Productivity  Advancement  Program).  The  formation 
of  a  consortium  to  integrate  the  defense  and  commercial  industrial  bases  is  of 
increasing  importance  and  is  sponsored  through  the  US  Government's 
Technology  Reinvestment  Project  for  Dual  Use  Technologies.  CRDA’s 
(Cooperative  Research  and  Development  Agreements)  can  be  used  to  facilitate 
cooperative  research  between  the  government  and  private  industry  while 
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protecting  the  interests  of  both.  The  research  of  the  future  will  need  to  include 
collaboration  between  industry,  government  and  academia  in  order  to  get  the 
highest  return  on  shrinking  research  dollars. 
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Computer  Simulation  Models  for  Evaluating  the 
Performance  and  Design  of  Tracked  and  Wheeled  Vehicles 

J.Y.  Wong1 

Summary 

A  series  of  computer  simulation  models  for  performance  and  design  evaluation  of 
tracked  vehicles  and  off-road  wheeled  vehicles  have  emerged  in  the  past  decade.  In  contrast 
with  empirical  models  developed  earlier,  they  are  based  on  detailed  studies  of  the  mechanics 
of  vehicle-terrain  interaction,  and  take  into  account  all  major  vehicle  design  features  and 
terrain  characteristics.  Thus,  they  provide  a  comprehensive  and  realistic  tool  for  the  vehicle 
engineer  to  optimize  vehicle  design  and  for  the  procurement  manager  to  evaluate  competing 
vehicle  candidates.  These  models  have  been  gaining  increasingly  wide  acceptance  in  industry 
and  governmental  agencies.  For  instance,  the  model  NTVPM  for  tracked  vehicles  with 
relatively  short  track  pitch  has  been  successfully  used  to  assist  vehicle  manufacturers  in  the 
development  of  a  new  generation  of  high-mobility  military  vehicles  and  governmental 
agencies  in  the  evaluation  of  vehicle  candidates  in  Europe,  North  America  and  Asia. 

Introduction 

In  the  past  two  decades,  a  variety  of  computer  simulation  models  (computer-aided 
methods)  for  evaluating  the  mobility  of  off-road  vehicles  have  emerged.  In  the  early  1970s, 
to  support  decision  making  processes  related  to  the  procurement  and  deployment  of  military 
vehicles,  an  empirical  model  known  as  the  U.S.  Army  Mobility  Model  (AMC-71)  was 
developed.  In  the  mid-1970s,  the  second  generation  of  this  model  called  AMM-75  was  made 
available  (Jurkat,  Nuttall,  and  Haley,  1975).  This  version  of  the  model  forms  the  basis  for  the 
subsequent  development  of  the  NATO  Reference  Mobility  Model  (NRMM).  NRMM  has  the 
capability,  among  others,  to  predict  the  tractive  performance  of  off-road  vehicles  over 
unprepared  terrain  using  empirical  relationships.  This  capability  is,  however,  limited  to  the 
prediction  of  vehicle  performance  over  two  types  of  terrain,  namely,  fine-  and  coarse-  grained 
soils. 

To  assist  the  development  and  design  engineer  to  optimize  the  design  of  off-road 
vehicles  and  the  procurement  manager  to  select  the  appropriate  vehicle  candidates  for  a  given 
mission  and  operating  environment,  a  series  of  computer  simulation  models  have  been 
developed  under  the  auspices  of  Vehicle  Systems  Development  Corporation  of  Canada,  since 
the  1980s.  In  contrast  with  empirical  models  developed  earlier,  these  models  are  based  on 
detailed  studies  of  the  physical  nature  of  vehicle-terrain  interaction  and  the  principles  of 
applied  mechanics.  They  take  into  account  all  major  design  features  of  the  vehicle  and  the 
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basic  characteristics  of  the  terrain. 

For  performance  and  design  evaluation  of  vehicle  with  flexible  tracks,  such  as  rubber- 
belt  tracks  or  link  tracks  with  relatively  short  track  pitch,  commonly  found  in  military  fighting 
and  logistics  vehicles,  a  model  known  as  NTVPM  has  been  developed  (Wong,  1986,  1989, 
1992a  and  b,  and  1993;  Wong  and  Preston-Thomas,  1986a  and  1988).  This  model  is  based 
on  the  assumption  that  this  type  of  track  can  be  idealized  as  a  flexible  belt.  It  has  been 
successfully  used  to  assist  manufacturers  in  the  development  of  new  military  vehicles  and 
governmental  agencies  in  the  evaluation  of  vehicle  candidates  in  Europe,  North  America  and 
Asia.  For  instance,  this  model  has  been  employed  to  assist  Hagglunds  Vehicle  AB  of  Sweden 
in  the  development  of  a  new  generation  of  high-mobility  fighting  vehicles  (CV  90),  in  the 
examination  of  the  approach  to  the  further  improvement  of  the  performance  of  the  all-terrain 
carrier  BV  206,  and  in  the  evaluation  of  competing  designs  for  a  proposed  main  battle  tank. 
Recently,  it  has  been  used  in  the  selection  of  an  optimum  configuration  for  a  new,  high- 
mobility  armoured  personnel  carrier  for  a  Spanish  vehicle  manufacturer  and  in  the  assessment 
of  the  effects  of  design  modifications  on  the  mobility  of  a  fighting  vehicle  over  tropical 
terrain  for  an  Asian  firm.  It  has  also  been  employed  in  the  evaluation  of  the  effects  of  design 
changes  on  the  cross-country  mobility  of  Canada’s  main  battle  tank.  Leopard  Cl,  for  the 
Canadian  Department  of  National  Defence  and  in  the  assessment  of  the  mobility  of  a  variety 
of  container  handling  equipment  used  by  the  U.S.  Marine  Corps. 

For  tracked  vehicles  with  rigid  links  having  relatively  long  track  pitch,  commonly  used 
in  agriculture  and  construction  industry,  a  model  known  as  RTVPM  has  recently  been 
developed  (Gao  and  Wong,  1993  and  1994;  Wong  and  Gao,  1994).  In  this  model,  the  track 
is  considered  to  be  a  system  of  rigid  links  connected  through  frictionless  pins.  The  basic 
features  of  this  model  have  been  verified  with  available  field  test  data. 

For  performance  and  design  evaluation  of  off-road  wheeled  vehicles,  a  model  known 
as  NWVPM  has  been  developed  (Wong  and  Preston-Thomas,  1986b).  It  can  be  employed  in 
the  evaluation  of  the  overall  performance  and  design  of  off-road  wheeled  vehicles,  as  well  as 
in  the  selection  of  tires  for  cross-country  operations.  It  has  been  used  in  the  assessment  of 
the  effects  of  different  types  of  tire  on  the  mobility  of  6  x  6  and  8x8  armoured  wheeled 
vehicles  for  the  Canadian  Department  of  National  Defence  and  in  the  evaluation  of  the 
mobility  of  container  handling  wheeled  vehicles  used  by  the  U.S.  Marine  Corps. 

In  this  paper,  the  basic  features  and  capabilities  of  these  computer  simulation  models 
will  be  presented.  Experimental  validations  of  these  models  with  field  test  data  will  be 
described.  The  applications  of  these  models  to  parametric  analysis  of  vehicle  performance 
and  to  the  optimization  of  vehicle  design  will  be  demonstrated. 

Computer  Simulation  Model  NTVPM  for  Vehicles  with  Flexible  Tracks 

For  high-speed  tracked  vehicles,  such  as  military  fighting  and  logistics  vehicles  and 
off-road  transport  vehicles,  rubber  tracks  or  link  tracks  with  relatively  short  track  pitch  are 
commonly  used.  This  kind  of  short-pitch  track  system  typically  has  a  ratio  of  roadwheel 
diameter  to  track  pitch  in  the  range  from  4  to  6,  a  ratio  of  roadwheel  spacing  to  track  pitch  in 


-  3  - 


the  range  of  4  to  7,  and  a  ratio  of  sprocket  pitch  diameter  to  track  pitch  usually  of  the  order 
of  4.  The  rubber-belt  track  and  the  short-pitch  link  track  will  be  referred  to  as  "flexible 
track"  in  this  paper  and  they  may  be  idealized  as  a  flexible  belt  in  the  analysis  of  track-terrain 
interaction. 

The  computer  simulation  model  NTVPM  has  been  developed  for  performance  and 
design  evaluation  of  tracked  vehicles  with  flexible  tracks,  under  the  auspices  of  Vehicle 
Systems  Development  Corporation.  The  model  is  intended  to  provide  the  vehicle  designer 
with  a  comprehensive  and  realistic  computer-aided  method  to  optimize  vehicle  design,  and  to 
provide  the  procurement  manager  with  a  reliable  method  to  evaluate  vehicle  candidates.  To 
meet  these  objectives,  the  latest  version  of  the  model,  known  as  NTVPM-86,  takes  into 
account  all  major  vehicle  design  parameters,  including  sprung  weight,  unsprung  weight, 
location  of  the  centre  of  gravity,  number  of  roadwheels,  location  of  roadwheels,  roadwheel 
dimensions  and  spacing,  locations  of  sprocket  and  idlers,  supporting  roller  arrangements,  track 
dimensions  and  geometry,  initial  track  tension,  belly  (hull)  shape,  and  angles  of  approach  and 
departure  of  the  track  system.  The  longitudinal  elasticity  of  rubber-belt  tracks  or  of  link 
tracks  with  rubber  bushings  are  taken  into  consideration.  The  track  longitudinal  elasticity 
affects  the  tension  distribution  in  the  track  and  as  a  result  influences  the  performance  of  the 
vehicle  to  a  certain  extent  over  marginal  terrain.  The  characteristics  of  the  independent 
suspension  of  the  roadwheels  are  fully  taken  into  account  in  the  model.  Torsion  bar 
suspensions,  hydro-pneumatic  suspensions  with  non-linear  load-deflection  characteristics,  and 
others  can  be  accommodated  in  the  model.  Suspensions  characteristics  have  a  significant 
effect  on  vehicle  mobility  over  soft  ground.  On  highly  compressible  terrain,  such  as  deep 
snow,  track  sinkage  may  be  greater  than  the  ground  clearance  of  the  vehicle.  If  this  occurs, 
the  belly  (hull)  of  the  vehicle  will  be  in  contact  with  the  terrain  surface  and  will  support  part 
of  the  vehicle  weight.  This  will  reduce  the  load  carried  by  the  tracks  and  will  adversely 
affect  the  traction  of  the  vehicle  over  terrain  that  exhibits  frictional  behaviour.  Furthermore, 
belly  contact  will  give  rise  to  an  additional  drag  component  -  the  belly  drag.  The  problem  of 
belly  contact  is  of  importance  to  vehicle  mobility  over  marginal  terrain,  and  the  characteristics 
of  belly-terrain  interaction  have  been  taken  into  consideration  in  the  model.  Terrain 
characteristics,  including  the  pressure-sinkage  relation,  shear  strength,  rubber-terrain  shearing 
(for  rubber  tracks  or  tracks  with  rubber  pads)  and  belly-terrain  shearing  characteristics,  and 
responses  to  repetitive  normal  and  shear  loadings,  are  taken  into  account  in  the  model. 

Basic  Approach  to  the  Development  of  the  Model 

In  developing  the  model,  the  track  is  assumed  to  be  a  flexible  belt  with  known 
longitudinal  elasticity.  The  track-roadwheel  system  used  in  the  analysis  is  schematically 
shown  in  Fig.  1.  When  a  tracked  vehicle  travels  over  a  deformable  terrain,  the  load  applied 
through  the  track  system  causes  the  terrain  to  deform.  The  track  segments  between 
roadwheels  take  up  load,  and  as  a  result  they  deflect  and  have  a  form  of  a  curve.  The  actual 
length  of  the  track  in  contact  with  the  terrain  between  the  front  and  rear  roadwheels  increases 
in  comparison  with  that  when  the  track  rests  on  a  firm  ground.  This  causes  a  reduction  in  the 
sag  of  the  top  run  of  the  track  and  a  change  in  track  tension  distribution.  It  should  also  be 
pointed  out  that  an  element  of  the  terrain  beneath  the  track  is  first  subject  to  the  load  applied 
by  the  leading  roadwheel.  When  the  leading  roadwheel  has  passed,  the  load  on  the  terrain 
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Figure  1.  A  flexible  track  model 
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element  is  relieved.  Load  is  reapplied  as  the  second  roadwheel  rolls  over  it.  A  terrain 
element  under  the  track  is  thus  subject  to  repetitive  loading.  The  loading-unloading-reloading 
cycle  continues  until  th'*  rear  roadwheel  of  the  track  system  has  passed  over  it.  To  predict 
the  normal  and  shear  stress  distributions  under  a  moving  tracked  vehicle,  the  pressure-sinkage 
relation,  shearing  characteristics,  and  responses  to  repetitive  normal  and  shear  loadings  of  the 
terrain  are  taken  into  consideration. 

Based  on  an  understanding  of  the  physical  nature  of  the  problem,  the  mechanics  of 
track-terrain  interaction  is  analyzed  in  detail.  A  set  of  equations  for  the  equilibrium  of  the 
forces  and  moments  acting  on  the  track  system  are  derived.  They  establish  the  relationship 
between  the  shape  of  the  deflected  track  in  contact  with  the  terrain  and  vehicle  design 
arameters  and  terrain  characteristics.  The  solutions  to  this  set  of  equations  define  the 
smkages  of  the  roadwheels,  the  inclination  of  the  vehicle,  the  track  tension  distribution,  and 
the  track  shape  in  contact  with  the  ground.  From  these,  the  normal  and  shear  stress 
distributions  on  the  track-terrain  interface,  and  the  track  motion  resistance,  belly  drag  (if  the 
vehicle  belly  is  in  contact  with  the  terrain),  thrust,  drawbar  pull,  and  tractive  efficiency  of  the 
vehicle  as  functions  of  track  slip  can  be  determined.  For  further  details  of  the  model,  please 
refer  to  the  references  (Wong,  1989  and  1993). 

Experimental  Validation 

The  model  can  be  used  to  predict  the  performance  of  single  unit  and  two-unit 
articulated  tracked  vehicles  over  unprepared  terrain.  Its  basic  features  have  been  validated 
with  field  test  data  obtaining  using  various  test  vehicles,  including  M113A1,  BV202  and 
BV206,  over  a  variety  of  unprepared  terrains,  including  mineral  terrain,  organic  terrain 
(muskeg)  and  snow-covered  terrain.  Figures  2-4  show  a  comparison  of  the  measured  and 
predicted  normal  pressure  distributions  under  the  track  pad  of  an  armoured  personnel  carrier 
M113A1  over  a  sandy  terrain,  a  muskeg,  and  a  snow-covered  terrain,  respectively.  A 
comparison  of  the  measured  and  predicted  drawbar  performance  of  an  M113A1  over  the  three 
types  of  terrain  is  shown  in  Figs.  5-7,  respectively.  Figure  8  shows  a  comparison  between 
the  measured  and  predicted  drawbar  performance  of  a  two-unit,  articulated  tracked  vehicle  BV 
206,  over  an  undisturbed  snow.  Reasonably  close  agreements  between  the  measured  and 
predicted  normal  pressure  distributions  and  drawbar  performance  obtained  using  NTVPM-86 
confirm  the  validity  of  the  basic  features  of  the  model. 

Applications  to  Parametric  Analysis  and  Design  Optimization 

NTVPM-86  can  be  employed  to  assess  the  effects  of  vehicle  design  parameters  on 
vehicle  mobility  and  the  influence  of  terrain  conditions  on  vehicle  performance.  The  model 
can  also  be  used  in  design  optimization  for  a  given  mission  and  operating  environment. 

Figure  9  shows  the  effects  of  the  number  of  roadwheels  and  the  initial  track  tension  on 
the  drawbar  pull  to  weight  ratio  (drawbar  pull  coefficient)  of  a  reference  vehicle  with  design 
parameters  similar  to  that  of  the  Ml  13 A1  over  a  deep  snow,  predicted  using  the  simulation 
model  (Wong  and  Preston-Thomas,  1986a).  It  can  be  seen  that  both  the  number  of 
roadwheels  and  the  initial  track  tension  have  significant  effects  on  vehicle  mobility  over  soft 
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Figure  5 .  Comparison  of  the  measured  and  predicted 
drawbar  performance  of  an  M113A1  on  a 
sandy  terrain 


Figure  6.  Comparison  of  the  measured  and  predicted 
drawbar  performance  of  an  M113A1  on  a 
muskeg 
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Figure  7 .  Comparison  of  the  measured  and  predicted 
drawbar  performance  of  an  M113A1 
on  a  snow 
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Figure  8.  Comparison  of  the  measured  and  predicted 
drawbar  performance  of  a  BV  206 
on  a  snow 
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Figure  10.  Effects  of  suspension  characteristics 
on  vehicle  performance  on  a  snow 


ground.  For  a  given  (or  existing)  vehicle,  its  mobility  over  marginal  terrain  can  be  greatly 
improved  by  increasing  the  initial  track  tension.  This  research  finding  obtained  using 
NTVPM-86  has  led  to  the  development  of  the  central  initial  track  tension  regulating  system 
controlled  by  the  driver.  Over  normal  terrain,  the  driver  can  set  the  track  tension  at  the 
regular  level.  However,  wnen  traversing  marginal  terrain  is  anticipated,  the  driver  can  readily 
increase  the  track  tension  to  an  appropriate  level  to  improve  vehicle  mobility.  The  central 
track  tension  regulating  system  is  analogous  to  the  central  tire  inflation  system  for  off-road 
wheeled  vehicles.  A  central  initial  track  tension  regulating  system  has  been  developed  and 
installed  on  a  new  generation  of  high-mobility  armoured  vehicles  currently  in  production  in 
Sweden. 


Figure  10  shows  the  effects  of  suspension  characteristics  on  the  mobility  over  deep 
snow  of  the  reference  vehicle  noted  above.  The  parameters  of  the  three  suspension 
configurations  examined  are  given  in  Table  1.  The  basic  difference  between  them  is  in  the 
settings  of  the  initial  torsion  arm  angles  under  no  load  conditions.  The  standard  configuration 
is  similar  to  that  of  the  Ml  13A1  with  the  initial  torsion  arm  angle  set  at  43°  for  all  roadwheel 
stations,  as  shown  in  Table  1.  For  suspension  configuration  S2,  the  initial  torsion  arm  angle 
is  set  in  a  decreasing  order  from  51.6°  at  the  front  (first)  roadwheel  station  to  34.4°  at  the 
rear  (fifth)  roadwheel  station,  while  maintaining  an  angle  of  43°  for  the  torsion  arm  at  the 
middle  (third)  roadwheel  station.  This  setting  results  in  a  nose-up  attitude  for  the  vehicle 
body.  In  deep  snow,  this  causes  the  load  supported  by  the  vehicle  belly  and  the  associated 
belly  drag  to  increase  and  vehicle  performance  to  decrease,  as  shown  in  Fig.  10.  For 
suspension  S3,  the  initial  torsion  arm  angle  is  set  in  an  increasing  order  from  34.4°  at  the 
front  (first)  roadwheel  station  to  51.6°  at  the  rear  (fifth)  roadwheel  station,  while  maintaining 
an  angle  of  43°  for  the  middle  (third)  roadwheel  station.  This  setting  results  in  a  nose-down 
attitude  for  the  vehicle  body.  In  a  deep  snow,  this  causes  a  reduction  in  the  belly  load  and 
belly  drag  and  hence  an  improvement  in  performance,  in  comparison  with  the  standard 
configuration  and  configuration  S2,  as  shown  in  Fig.  10. 


Table  1.  Torsion  Arm  Settings  for  the  Standard  Suspension  and  Suspensions  S2  and  S3 


Roadwheel 

Station 

Initial  Torsion  Arm  Angles  Under  No  Load 
(below  the  horizontal),  degrees 

Suspension  Configuration 

Standard 

S2 

S3 

1 

43 

51.5 

34.4 

2 

43 

47.3 

38.7 

3 

43 

43 

43 

4 

43 

38.7 

47.3 

5 

43 

34.4 

51.6 

The  model  NTVPM-86  can  also  be  used  for  design  optimization  of  tracked  vehicles 
with  flexible  tracks.  Figure  1 1  shows  the  drawbar  performance  over  deep  snow  of  four 
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vehicle  configurations.  Configurations  A  and  B,  the  standard  configuration  with  parameters 
similar  to  that  of  the  Ml  13A1,  and  a  vehicle  configuration  similar  to  the  standard  one  but 
with  an  initial  track  tension  to  vehicle  weight  ratio  of  0.3.  Configuration  A  has  the 
suspension  configuration  S3  described  in  Table  1,  an  initial  track  tension  to  vehicle  weight 
ratio  of  0.25,  a  track  width  of  75  cm  and  a  ground  clearance  of  52  cm.  Configuration  B  has 
the  suspension  configuration  S3,  an  initial  track  tension  to  vehicle  weight  ratio  of  0.3,  a  track 
width  of  100  cm  and  a  ground  clearance  of  57  cm.  It  is  shown  that  Configurations  A  and  B 
exhibit  superior  tractive  performance  in  deep  snow  over  the  standard  configuration.  This 
indicates  that  NTVPM-86  can  be  an  extremely  useful  tool  for  the  design  engineer  to  evaluate 
competing  vehicle  designs  and  to  select  the  optimum  configuration  for  given  operating 
requirements. 

Computer  Simulation  Model  RTVPM  for  Vehicles  with  Rigid  Link  Tracks 

For  low-speed  tracked  vehicles,  such  as  those  used  in  agriculture  and  construction 
industry,  rigid  link  tracks  with  relatively  long  track  pitch  are  commonly  used.  This  type  of 
track  system  has  a  ratio  of  roadwheel  diameter  to  track  pitch  as  low  as  1.2  and  a  ratio  of 
roadwheel  spacing  to  track  pitch  typically  1.5. 

The  computer  simulation  model  RTVPM  has  been  developed  for  performance  and 
design  evaluation  of  tracked  vehicles  wth  rigid  link  tracks.  This  model  takes  into  account  all 
major  design  parameters  of  the  vehicle,  including  vehicle  weight,  location  of  the  centre  of 
gravity,  number  of  roadwheels,  location  of  roadwheels,  roadwheel  dimensions  and  spacing, 
locations  of  sprocket  and  idlers,  supporting  roller  arrangements,  track  dimensions  and 
geometry,  initial  track  tension,  and  drawbar  hitch  location.  As  the  track  links  are  considered 
to  be  rigid,  the  track  is  assumed  to  be  inextensible.  For  most  low-speed  tracked  vehicles,  the 
roadwheels  are  not  sprung,  and  hence  considered  to  be  rigidly  connected  to  the  track  frame. 
Terrain  parameters  used  in  this  model  are  the  same  as  those  used  in  NTVPM. 

Basic  Approach  to  the  Development  of  the  Model 

The  model  RTVPM  treats  the  track  as  a  system  of  rigid  links  connected  with 
frictionless  pin,  as  shown  in  Fig.  12.  As  noted  previously,  the  roadwheels,  supporting  rollers, 
and  sprocket  are  assumed  to  be  rigidly  attached  to  the  vehicle  frame.  The  centre  of  the  front 
idler  is,  however,  assumed  to  be  mounted  on  a  pre-compressed  spring. 

In  the  analysis,  the  track  system  is  divided  into  four  sections:  the  upper  run  of  the 
track  supported  by  rollers;  the  lower  run  of  the  track  in  contact  with  the  roadwheels  and  the 
terrain;  the  section  in  contact  with  the  idler;  and  the  section  in  contact  with  the  sprocket.  By 
considering  the  equilibrium  of  various  sections  of  the  track  system,  the  interaction  between 
the  lower  run  of  the  track  and  the  terrain,  and  the  compatibility  conditions  for  various  track 
sections,  a  set  of  equations  can  be  formulated.  The  solutions  to  this  set  of  equations 
determine  the  sinkage  and  inclination  of  the  track  system,  the  normal  and  shear  stress 
distributions  on  the  track-terrain  interface,  and  the  track  motion  resistance,  thrust,  drawbar 
pull,  and  tractive  efficiency  of  the  vehicle  as  functions  of  track  slip.  Figure  13  shows  the 
predicted  normal  pressure  distribution  under  a  tracked  vehicle  with  seven  roadwheels  on  a 
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Figure  13.  The  normal  pressure  distribution  under 
a  tractor  with  rigid  links  predicted 
using  RTVPM  on  a  clayed  soil 


Figure  14.  Comparison  of  the  measured  and  predicted 
drawbar  performance  of  a  tractor  on  a 
sandy  loam 
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clayey  soil  (Wong  and  Gao,  1994).  The  detailed  description  of  the  model  may  be  found  in 
the  references  (Gao  and  Wong,  1993  and  1994). 

Experimental  Validation 

The  basic  features  of  RTVPM  have  been  validated  with  available  field  test  data. 
Figures  14  and  15  show  a  comparison  of  the  measured  and  predicted  drawbar  pull  coefficient 
and  tractive  efficiency  as  functions  of  track  slip,  respectively,  for  a  heavy  tracked  vehicle 
used  in  construction  industry.  The  vehicle  has  a  total  weight  of  329  kN.  It  has  eight 
roadwheels  of  diameter  26  cm  on  each  of  the  two  tracks  and  the  average  spacing  between 
roadwheels  is  34  cm.  The  track  pitch  is  21.6  cm  and  the  track  width  is  50.8  cm.  The  terrain 
is  a  dry,  disked  sandy  loam,  with  an  angle  of  shearing  resistance  of  40.1°  and  a  cohesion  of 
0.55  kPa.  The  measured  data  shown  in  figures  are  provided  by  Caterpillar  Inc.,  Peoria, 
Illinois,  U.S.A. 

It  can  be  seen  that  the  tractive  performance  of  the  vehicle  predicted  using  RTVPM  is 
very  close  to  the  measured  one.  This  suggests  that  the  model  is  capable  of  providing  realistic 
predictions  of  vehicle  performance  in  the  field.  It  would  be  desirable,  however,  to  further 
validate  the  model  over  a  wider  range  of  terrain  conditions. 

Applications  to  Parametric  Analysis  and  Design  Optimization 

The  applications  of  RTVPM  to  design  evaluation  and  parametric  study  of  vehicles 
with  rigid  link  tracks  will  be  demonstrated  through  examples. 

Figure  16  shows  the  effects  of  the  track  pitch  on  the  drawbar  performance  on  a  clayey 
soil  of  a  reference  vehicle  with  a  total  weight  of  372.4  kN,  predicted  using  RTVPM  (Wong 
and  Gao,  1994).  The  vehicle  has  seven  roadwheels  on  each  track,  a  track  pitch  of  21.6  cm, 
an  initial  track  tension  of  22.30  kN,  and  a  centre  of  gravity  at  the  mid-point  of  the  track 
contact  length.  It  can  be  seen  that  within  the  range  studied,  the  longer  the  track  pitch,  the 
higher  the  tractive  performance  will  be.  This  is  primarily  due  to  the  fact  that  with  a  longer 
track  pitch,  the  normal  pressure  distribution  under  the  track  becomes  more  favourable.  It 
should  be  noted,  however,  that  with  a  longer  track  pitch,  the  fluctuation  in  speed  and  the 
vibrations  of  the  track  system  may  increase.  Consequently,  there  is  a  practical  limit  to  the 
track  pitch  for  a  given  vehicle  configuration. 

Figure  17  shows  the  effects  of  the  number  of  roadwheels  on  the  drawbar  performance 
of  the  reference  vehicle  on  the  clayey  soil.  It  can  be  seen  that  increasing  the  number  of 
roadwheels  enhances  the  tractive  performance  of  the  vehicle.  The  improvement  in  tractive 
performance  with  a  larger  number  of  roadwheels  is  due  to  a  more  uniform  normal  pressure 
distribution. 

The  model  can  also  be  used  for  design  optimization  of  tracked  vehicles  with  rigid  link 
'racks.  Figure  18  shows  a  comparison  of  the  drawbar  performance  of  the  reference  vehicle 
and  an  optimized  configuration  (Configuration  A)  on  the  clayey  soil.  Configuration  A  has 
nine  roadwheels  on  each  track,  a  track  pitch  of  23  cm,  an  initial  track  tension  of  89.20  kN, 
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Figure  16 .  Effects  of  track  pitch  on  the  performance 
of  a  tractor  on  a  clayey  soil 
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Figure  17.  Effects  of  the  number  of  roadwheels  on  the 
performance  of  a  tractor  on  a  clayey  soil 
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and  a  centre  of  gravity  at  40  cm  ahead  of  the  mid-point  of  the  track  contact  length.  It  shows 
that  RTVPM  can  be  a  useful  tool  for  the  vehicle  designer  to  select  the  optimum  vehicle 
configuration  and  design  parameters. 

Computer  Simulation  Model  NWVPM  for  Off-Road  Wheeled  Vehicles 

NWVPM  has  been  developed  for  the  evaluation  of  the  overall  performance  and  design 
of  off-road  wheeled  vehicles  over  unprepared  terrain,  as  well  as  for  the  selection  of  tires  for 
cross-country  operations.  The  model  takes  into  account  all  major  design  parameters  of  the 
vehicle  as  well  as  the  tire.  The  vehicle  design  parameters  considered  include:  vehicle  weight, 
axle  load,  axle  spacing,  location  of  the  centre  of  gravity,  axle  suspension  stiffness,  function  of 
the  axle  (driven  or  non-driven),  axle  clearance,  track  of  the  axle,  belly  (hull)  shape,  and 
drawbar  hitch  location.  The  tire  parameters  considered  include:  outside  diameter,  tread 
width,  section  height,  lug  area/carcass  area,  lug  height,  lug  width,  inflation  pressure,  average 
ground  contact  pressure,  and  tire  construction  (radial  or  bias).  Terrain  parameters  used  in  the 
model  are  the  same  as  those  used  in  NTVPM  and  RTVPM. 

Basic  Approach  to  the  Development  of  the  Model 

The  model  NWVPM  consists  of  two  sub-models,  one  is  the  tire  sub-model  and  the 
other  is  the  vehicle  sub-model. 

The  tire  sub-model  used  is  that  developed  by  Wong  (Wong,  1989  and  1993).  Based 
on  the  dimensions  of  the  tire,  the  inflation  pressure  and  carcass  stiffness  (or  alternatively  the 
average  contact  pressure  on  a  hard  surface),  the  normal  load,  and  terrain  characteristics,  the 
operating  mode  of  the  tire  ("rigid"  or  "elastic")  is  first  predicted.  Based  on  a  detailed  analysis 
of  the  mechanics  of  tire-terrain  interaction,  a  set  of  equations  for  the  equilibrium  of  the  tire 
can  then  be  established.  The  solutions  to  this  set  of  equations  determine  the  normal  and  shear 
stress  distributions  on  the  tire-terrain  interface,  the  motion  resistance,  (including  the  internal 
resistance  of  the  tire),  thrust,  and  sinkage  of  the  tire.  A  schematic  of  the  tire  sub-model  used 
is  shown  in  Fig.  19. 

The  tire  sub-model  is  incorporated  into  the  vehicle  sub-model  to  provide  a  complete 
framework  for  performance  and  design  evaluation  of  off-road  wheeled  vehicles.  The  vehicle 
sub-model  takes  into  account  the  dynamic  inter-axle  load  transfer  and  the  suspension  stiffness 
of  the  axles.  Any  number  of  axles  can  be  accommodated.  When  the  track  of  the  front 
(preceding)  axle  is  the  same  as  that  of  the  rear  (following)  axle,  the  tires  on  the  rear  axle  run 
in  the  ruts  formed  by  the  tires  on  the  front  axle.  Terrain  properties  in  the  rut  will  be  different 
from  those  in  the  virgin  state.  To  take  into  account  this  "multipass"  effect,  the  responses  of 
the  terrain  to  repetitive  normal  and  shear  loadings  are  taken  into  consideration  in  the  model. 

In  addition,  both  single  and  dual  tires  can  be  accommodated.  The  output  of  the  model 
includes  the  load,  sinkage,  motion  resistance  and  thrust  of  the  axles,  and  the  drawbar  pull  and 
tractive  efficiency  of  the  vehicle  as  functions  of  wheel  slip. 

The  basic  features  of  the  model  have  been  verified  with  available  field  test  data. 
Figures  20  and  21  show  a  comparison  of  the  measured  and  predicted  drawbar  performance  of 
a  tractor  obtained  using  NWVPM  on  a  plowed  and  stubble  field,  respectively. 
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Figure  19.  A  tire  model 
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Applications  to  Parametric  Analysis  and  Design  Optimization 

The  applications  of  NWVPM  to  parametric  analysis  of  the  performance  and  design  of 
off-road  wheeled  vehicles  and  to  the  selection  of  tires  for  a  given  operating  environment  will 
be  illustrated  through  examples. 

Figure  22  shows  the  effects  of  tire  design  and  inflation  pressures  of  the  front  and  rear 
tires  on  the  tractive  performance  of  a  two-axle  vehicle  on  a  medium  soil,  predicted  using 
NWVPM.  The  first  and  second  numbers  in  the  inflation  pressure  combinations  shown  in  the 
figure  represent  the  inflation  pressure  of  the  front  tires  and  that  of  the  rear  tires,  respectively. 
The  effects  of  the  static  load  distribution  between  the  front  and  rear  axles  on  the  drawbar 
performance  of  the  two-axle  vehicle  are  shown  in  Fig.  23.  It  indicates  that  because  of  the 
"multipass"  effect,  a  lighter  static  load  on  the  front  and  a  heavier  static  load  on  the  rear  will 
give  improved  tractive  performance  on  the  medium  soil  (Wong,  1989). 

Concluding  Remarks 

Computer  simulation  models  based  on  empirical  relations  have  played  a  useful  role  in 
the  past.  However,  with  an  improved  understanding  of  the  physical  nature  of  vehicle-terrain 
interaction  and  of  terrain  response  to  vehicular  loading,  a  new  generation  of  computer 
simulation  models  has  emerged  over  the  past  decade.  They  are  based  on  detailed  analyses  of 
the  mechanics  of  vehicle-terrain  interaction  and  take  into  account  all  major  vehicle  design 
features  and  terrain  characteristics.  These  comprehensive  and  realistic  computer  simulation 
models  have  played  and  will  continue  to  play  an  increasingly  important  role  in  the  future 
development  of  off-road  vehicles.  For  instance,  the  computer  simulation  model  NTVPM  for 
vehicles  with  flexible  tracks  have  been  successfully  used  to  assist  vehicle  manufacturers  in  the 
development  of  new  products  and  governmental  agencies  in  the  evaluation  of  vehicle 
candidates  in  Europe,  North  America  and  Asia. 

To  further  develop  computer  simulation  models  (computer-aided  methods)  for 
performance  and  design  evaluation  of  off-road  vehicles,  the  following  guidelines  are 
suggested. 

a)  It  should  be  clear  that  the  objective  of  a  model  is  to  provide  a  framework  that  will 
enable  the  engineering  practitioner  to  reaho  .ically  evaluate  the  performance  or  design 
of  off-road  vehicles  under  a  variety  of  operating  environments.  The  development  of 
the  model  is  an  engineering  endeavour  and  not  an  academic  exercise. 

b)  The  model  should  address  the  needs  of  the  vehicle  user,  designer,  or  procurement 
manager  and  not  that  of  the  theoretician. 

c)  The  model  should  be  developed  and  implemented  in  such  a  manner  that  will  be 
conducive  to  practical  results  and  should  appeal  to  a  wide  spectrum  of  engineering 
practitioners. 
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Figure  23.  Effects  of  static  weight  distribution 
on  the  performance  of  a  two-axle 
wheeled  vehicle  on  a  medium  soil 


d)  In  the  development  of  the  model,  including  the  characterization  of  terrain  behaviour,  a 
pragmatic  engineering  approach  should  be  followed.  It  should  not  be  unduly  dwelling 
on  theoretical  niceties. 

e)  An  off-road  vehicle  is  a  complex  mechanical  system.  In  the  development  of  the 
model,  emphasis  should  be  placed  on  an  adequate  representation  of  the  vehicle,  so  that 
meaningful  results  can  be  obtained  to  guide  its  development,  design  and  procurement. 

f)  While  an  understanding  of  the  mechanical  behaviour  of  the  terrain  (soil)  is  essential, 
the  development  of  the  model  is  not  an  exercise  in  soil  mechanics. 

g)  The  success  or  failure  of  a  model  is  eventually  judged  by  the  market  place  or  by  the 
engineering  practitioner,  and  not  by  the  theoretician  or  the  bureaucrat. 
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SUMMARY: 

Soil-tire/Soil-track  interaction  is  of  particular  interest  to  researchers  involved  in  off-road 
mobility  and  traction  research.  This  includes  scientists  and  engineers  involved  in  research 
in  the  field  of  agriculture,  construction,  forestry,  military,  and  mining.  In  agriculture  and 
forestry  soil  compaction  caused  by  traction  devices  is  also  a  serious  concern.  A  sound 
mathematical  model  is  a  pre-requisite  to  obtain  a  clear  understanding  of  the  soil-tire/soil- 
track  interaction  process.  A  key  ingredient  for  any  such  model  is  a  constitutive  relationship 
which  describes  the  stress-strain  behavior  of  soil.  Any  suitable  constitutive  model  requires 
soil  physical  properties  which  describe  the  elastic  behavior  of  soil,  onset  of  yield  and 
subsequent  plastic  flow,  material  hardening  or  softening  rules  etc.  Since  in-situ  soils 
seldom  behave  like  remolded  laboratory  soils  or  disturbed  field  samples,  it  is  important  to 
"identify"  or  "calibrate"  the  engineering  properties  of  field  soil  by  means  of  in-situ  tests. 
The  technique  of  obtaining  material  parameters  based  on  actual  system  response  is  known 
as  "back  analysis",  "inverse  solution",  "identification",  or  "calibration  procedure".  For 
complex  problems  such  as  soil-traction  device  interaction  where  closed  form  analytical 
solutions  do  not  exist  a  numerical  technique  such  as  a  finite  element  technique  is 
commomnly  used  to  solve  underlying  system  differential  equation.  For  such  cases  the 
back  analysis  procedure  can  take  one  of  the  two  forms:  (1)  inverse  method,  and  (2)  direct 
method.  This  paper  addresses  the  advantages  and  disadvantages  of  such  techniques,  and 
discusses  a  new  technique  which  overcomes  some  of  their  limitations.  This  new  technique 
consists  of  developing  a  so  called  "response  surface"  in  the  parameter  space  and  then  using 
this  pie-determined  surface  to  "identify"  engineering  properties  of  the  material  based  on 
in-situ  tests.  Two  case  studies  -  (1)  a  two  parameter  hypo-elastic  model  for  soil,  and  (2)  a 
complex  five  parameter  model  for  soil  which  includes  nonlinear  material  behavior  in  elastic 
range,  yield  based  on  Drucker-Prager  yield  criteria  and  associated  plastic  flow  upon  yield 
are  presented  to  illustrate  the  methodology. 


INTRODUCTION  AND  REVIEW  OF  LITERATURE: 

One  of  the  challenges  in  the  design  of  an  off-road  vehicle  is  to  equip  it  with  a  traction 
devicef  tire  or  track)  which  can  develop  high  traction  efficiently(  i.e.  optimum  tractive 
efficiency)  while  deterring  soil  compaction.  Even  an  increase  of  one  percentage  point  in 
the  tractive  efficiency  leads  to  an  annual  savings  of  over  100  million  liters(  about  25  million 
gallons)  of  fuel  in  U.  S.  alonefl].  On  the  other  hand,  soil  compaction  has  been  recognized 
as  a  worldwide  problem  with  serious  implications  on  agricultural  sustainability [2]. 
Although,  certain  amount  of  soil  compaction  may  even  be  desirable  for  some  crops  under 
certain  environmental  conditions  (optimum  soil  compaction),  excessive  soil  compaction 
can  lead  to  diminished  soil  porosity,  reduced  water  infiltration,  increased  resistance  to  root 
penetration,  increased  tillage  energy  requirements,  decreased  biological  activity,  and  a 
reduction  in  crop  yield[3  - 14].  A  necessary  pre-requisite  for  the  successful  design  of  a 
traction  device  is  a  sound  mathematical  model  for  the  soil-traction  interaction  process.  This 
interaction  is  an  extremely  complex,  dynamic  process.  A  key  ingredient  of  such  a  model  is 
a  constitutive  relationship  which  describes  the  stress-strain  behavior  for  soil.  Schafer  et  al. 
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[15]  stated  that  an  accurate  description  of  soil  constitutive  relationship  is  necessary  for  the 
integrity  and  robustness  of  the  model.  Soil  is  perhaps  one  of  the  most  complex  material 
from  engineering  point  of  view[16]. 

Numerous  constitutive  models  are  currently  available  for  soils.  Among  these  are  the 
elasticity  models,  higher  order  nonlinear  elasticity  models,  hypo-elasticity  models, 
plasticity  models  and  visco-plasticity  models.  Desai  [16],  Desai  and  Siriwardane  [17]  and 
Chen  and  Baladi[18]  have  discussed  these  models  and  their  applicability  to  a  specific 
loading  situation  in  detail.  Piece-wise  linear  elastic  models  (hyperbola,  parabola,  splines 
and  Ramberg-Osgood  formulas)  tend  to  be  good  for  a  specific  loading  case  but  are  poor  to 
simulate  general  loading  conditions.  Uigher  order  nonlinear  elasticity  models  tend  to 
include  too  many  parameters  and  have  limited  appeal.  Hypo-elasticity  models  appear  to 
show  some  promise.  Plasticity  models  which  utilize  Von  Mises,  Mohr-Coulomb  and 
Drucker-Prager  failure  criteria  have  been  widely  used.  To  include  volume  changes  due  to 
shear  in  geological  materials  and  also  to  account  for  strain  hardening  or  softening  behavior 
critical  state  models  have  been  developed.  CAM  and  CAP  models  account  for  growth  of  the 
yield  surface  and  have  become  increasingly  popular  in  civil  engineering.  Applicability  of 
critical  state  models  to  unsaturated  agricultural  soils  has  been  a  much  debated  issue. 
Hettiaratchi  and  O'  Callaghan[19],  Hettiaratchi[20]  and  Kirby[21]  have  found  that  critical 
state  concept  is  applicable  to  unsaturated  soils  both  qualitatively  and  quantitatively  except 
that  the  critical  state  parameters  depend  on  the  soil  moisture  content  They  found  that  it  is 
reasonable  to  use  total  stress  in  the  model(i.e.  soil  moisture  tension  can  be  ignored).  Bailey 
et  al.[22]  and  Bailey  and  Johnson[  23]  developed  a  constitutive  model  for  agricultural  soil 
that  relates  volumetric  strain  to  octahedral  normal  and  shear  stress.  This  model  predicts 
volumetric  strain  of  soil  samples  accurately  at  limiting  values  of  stresses(i.e.  zero  and  very 
large  applied  stress).  Raper  and  Erbach[24]  and  Raper  et  al.[25]  have  used  this  constitutive 
equation  to  compute  tangent  moduli  in  a  finite  element  program  to  predict  soil  compaction. 

All  the  aforementioned  constitutive  models  require  material  parameters.  These  material 
properties  describe  the  elastic  behavior  of  soil,  onset  of  yield  and  subsequent  plastic  flow, 
material  hardening  or  softening  rules  etc.  Typically  these  parameters  are  determined  using 
laboratory  tests.  Sometimes  remolded  soils  are  employed  in  the  laboratory  tests  which 
may  not  behave  like  field  soil.  Use  of  soil  properties  obtained  from  remolded  samples  can 
often  lead  to  predictions  which  are  unrealistic  and  of  little  value  to  engineers  interested  in 
improving  tire  design.  Even  if  field  samples  are  obtained,  one  of  the  main  problem  with 
the  soil  material  is  that  these  samples  undergo  disturbances  during  excavation  and  testing, 
and  may  not  behave  like  in-situ  soil  under  actual  loading  conditions  in  the  field.  Use  of 
c  ^  penetrometer,  grouser  plate,  and  sinkage  plate  often  yield  some  composite  soil 
parameters  which  depend  on  the  geometry  of  the  test  device  and  loading  conditions.  These 
composite  soil  parameters  are  of  little  use  in  subsequent  model  studies  based  on 
constitutive  relationship.  It  is  preferable  to  determine  the  soil  material  parameters  based  on 
undisturbed  in-situ  tests.  The  technique  of  obtaining  material  parameters  based  on  actual 
system  response  is  called  "back  analysis",  "inverse  solution",  "identification",  or  " 
calibration  procedure".  The  process  of  "calibrating"  actual  field  response  to  model 
behavior  is  expected  to  "identify"  the  material  parameters  which  can  accurately  predict 
system  response  in  subsequent  analysis  which  utilize  the  same  constitutive  model. 

The  back  analysis  technique  has  been  successfully  used  in  Geomechanics  in  studying 
tunneling  problems  in  rocks  and  in  investigating  settlement  problems[26-43].  If  a  closed 
form  solution  exists  for  the  underlying  differential  equation  describing  the  physical 
problem,  then  back  analysis  to  obtain  the  material  parameters  involves  optimizing  the 
difference  between  the  analytical  and  experimental  responses.  However,  most  real  life 
problems  in  geomechanics  are  geometrically  and/or  materially  nonlinear,  and  an  analytical 
solution  may  not  exist.  In  such  cases  a  numerical  procedure  such  as  a  finite  element 
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method[FEM]  may  be  used  to  obtain  solutions  to  the  governing  differential  equation. 
When  finite  element  analysis  is  used,  back  analysis  may  take  one  of  the  two  forms  -  1) 
inverse  method,  and  2)  direct  method. 

In  the  inverse  method  nodal  values  of  displacements  and  stresses  obtained  by  a  FEM 
technique  are  used  as  known  boundary  conditions  and  the  unknown  displacements  and 
stresses  are  eliminated  from  the  global  matrix  equation  by  reduction[41].  A  brief  discussion 
of  the  method  is  as  follows: 

Let  the  FEM  result  in  the  following  matrix  equation: 

(K]{u)  =  {F}  (1) 

where  K  is  the  global  stiffness  matrix,  u  is  the  nodal  displacement  vector  and  F  is  the 
global  forcing  vector.  Let  us  partition  the  global  stiffness  matrix  by  collecting  all  nodes  at 
which  nodal  values  are  measured  in  the  field  as  follows: 


‘K„ 

Ki2‘ 

“I 

►  —  < 

Fi’ 

_K2i 

k22. 

.U2. 

F2  J 

(2) 


where  u*i  is  a  vector  containing  measured  nodal  values  and  U2  is  a  vector  containing 
unknown  nodal  values.  Fi  and  F2  are  known  nodal  force  vectors,  and  Kij ’s  [  i=  1,2;  j= 
1,2  ]  are  partitioned  global  stiffness  matrix  elements.  Note  Ky ’s  are  functions  of  unknown 
material  parameter  vector,  p'.  Eliminating  U2  out  through  reduction,  we  get 

[K*]{u*>  =  {F*}  (3) 

where 

[K'l  =  [K„  +  K,jK- 'Kj,] 

{□•}  =  {□;} 

{F*}  =  {F,  -  K.jKj-j'Fj} 


In  equation  (3)  only  unknowns  are  pi  s  contained  in  the  elements  of  matrix  [K*].  An 
iterative  scheme  or  a  least  square  optimization  scheme  can  be  used  to  solve  equation  for 
unknown  material  parameters.  This  inverse  technique  is  quite  sensitive  to  experimental 
error  and  may  not  converge  at  all  in  some  cases[35, 40,42],  The  direct  approach  results  in 
more  accurate  parameter  values.  In  the  direct  method,  nodal  values  of  the  response  are 
computed  using  a  finite  element  method  for  a  set  of  assumed  parameter  values.  These 
responses  are  a  function  of  assumed  parameter  vector  (p),  say  u(£).  The  actual  values  of 
response  at  the  same  nodes  can  be  obtained  by  field  or  in-situ  tests.  If  u*  is  corresponding 
observed  response  to  u(ji)  then  ej=(uY  u(p)j)  is  a  measure  of  error  in  the  1th  value.  A 
suitable  objective  function  such  as  <j>=L  ei2  can  be  optimized  using  a  nonlinear  optimization 
technique[35,40,42].  The  direct  search  methods  such  as  simplex  method  or  its  modification 
such  as  Rosenbrock's  version  or  gradient  based  methods  such  as  conjugate  gradient 
method  or  quasi-Newton  method  can  be  successfully  used  depending  on  the 
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application[26,32,43].  Nodal  displacement  values  are  usually  better  than  stress  values  in 
parameter  identification[32,35,36].  Moreover,  it  is  preferable  to  map  all  the  parameters  to 
same  range  through  scaling[32].  Even  in  the  case  of  simple  linear  elastic  constitutive 

model,  the  objective  function,  <J>  will  be  a  nonlinear  function  of  material  parameter 

vector^  Because  of  this  situation,  the  objective  function,  4>  may  have  several  local 
miniraas[43].  Therefore,  the  optimal  solution  may  be  sensitive  to  initial  guess  values. 
Sometimes  different  combination  of  two  or  more  parameters  may  lead  to  same 
response[non-unique  solution][43].  More  than  one  type  of  test  or  tests  using  different 
geometry  and /  or  loads  may  be  helpful  in  such  cases.  Bayesian  approach  and  Kalman 
filtering  have  been  found  to  be  helpful  in  improving  the  accuracy  of  results  in  the  presence 
of  experimental  errors[27,33,40j.  The  direct  method  can  be  computationally  very 
expensive  since  at  each  iteration  a  new  FEM  analysis  with  updated  parameter  vector(p) 
needs  to  be  carried  out[42]. 


Rubinstein,  Upadhyaya,  and  Sime[44]  proposed  a  new  methodology  which  utilized 
orthogonal  regression  technique  to  develop  a  response  surface  in  the  parameter  space 
based  on  an  analytical  or  numerical  (such  as  a  finite  element  analysis)  solution  to  the 
system  differential  equation.  This  response  surface  was  used  in  the  optimization  step. 
Their  methodology  consists  of  following  steps: 


1.  A  response  surface  is  built  using  an  orthogonal  regression  technique  based  on  an 
analytical  or  numerical  solution  to  the  governing  differential  equation  of  the  system. 
The  response  surface  will  be  a  function  of  unknown  material  parameters. 

2.  This  response  surface  is  updated  using  higher  order  corrections  so  that  the  response 
surface  behavior  is  close  to  the  real  surface  behavior  everywhere  in  the  parameter  space. 
This  response  surface  will  be  used  to  predict  the  response  corresponding  to  the 
experimental  values(i.e.  at  the  same  load  and  nodal  point). 

3 .  Experimental  results  are  transformed  such  that  the  real  surface  and  the  response  surface 
will  have  one-to-one  correspondence  everywhere  in  the  region. 

4.  Experimental  results  are  optimized  against  the  response  surface  predictions  to  obtain 
material  properties  of  the  test  material. 


The  proposed  technique  is  particularly  useful  in  dealing  with  complex  problems  which 
require  numerical  solution  such  as  a  FEM  solution  to  the  underlying  system  differential 
equation.  The  main  advantage  of  this  technique  is  that  once  the  response  surface  is 
created  using  an  FEM  analysis,  there  is  no  need  to  go  back  to  the  FEM  analysis.  In  the 
classical  direct  or  indirect  approach,  hundreds  or  even  thousands  of  time  consuming  and 
expensive  FEM  evaluations  are  necessary  to  determine  material  parameters  through 
optimization  technique.  In  this  methodology  during  the  optimization  technique  only  the 
response  surface  is  used  to  estimate  u(£).  This  approach  is  expected  to  make  this 
technique  computationally  very  efficient  These  in-situ  soil  properties  can  be  used  in 
subsequent  model  studies  based  on  constitutive  relationships  which  utilize  these  soil 
parameters.  In  fact,  the  methodology  is  quite  general  and  can  be  used  in  other  fields  to 
estimate  constitutive  equation  parameters  based  on  in-situ  measurements. 


5 


MATHEMATICAL  MODELING 
Response  Surface  Development : 


Let  us  consider  a  general  material  constitutive  model  for  soil  (or  any  other  material) 
consisting  of  m  parameters:  PitP2.P3.~- .pm-  For  example,  if  we  select  a  nonlinear 
constitutive  model  with  extended  Drucker-Prager  yield  criteria  and  associated  flow  rule, 
then  six  parameters  will  be  involved[45,46].  These  parameters  are  p ^logarithmic  bulk 

rao(  k;  p2=  Poisson's  ratio,  t>;  p3=yield  surface  shape  factor(i.e.  related  to  the  third 

invariant  of  stress),  K;  p4=cohesion,  c;  p5=  internal  angle  of  friction,  <|>;  p6=initial  void 
ratio,  e.  The  last  parameter,  e  is  really  related  to  initial  stress  condition.  The  response  of  a 
system  to  applied  load  depends  on  its  geometry,  material  properties  and  the  load  itself.  If 
the  applied  load  and  the  geometry  are  fixed(i.e.  for  a  given  geometry  and  loading),  the 
system  response  is  a  function  of  material  constants  used  in  the  constitutive  equation. 

There  is  a  function  0=0(p1,p2,p3,...,pm)  which  represents  the  system  response  as  the 
material  properties  used  in  the  constitutive  equation  are  changed.  In  most  real  situations 
the  differential  equation  describing  the  response  is  nonlinear,  this  function  is  seldom 
known  explicitly.  One  of  the  goals  of  this  study  is  to  find  an  approximate  representation 

for  this  real  response,  <I>.  This  approximation  to  the  real  response  is  termed  the  response 
surface,  F  in  this  study.  One  convenient  way  of  determining  the  response  surface  F  is  to 
determine  the  variation  of  F  as  one  of  the  material  parameter,  pj  is  changed  while  all  other 
parameters  are  held  constant  Let  this  response  function  for  the  single  variable  pj  be  fj(pi). 

If  we  repeat  this  process  for  each  of  the  m  material  parameters(i.e.  for  i=l,2 . m),  then 

one  easy  way  of  obtaining  the  response  surface  is  simply  to  multiply  these  component 
equations,  f(pj),  i.e. 


<D  =  F  =  Cf1(p1)f2(p2)f3(P3) . fm(Pm)  ( 4) 

where 

F  =  response  surface 

f ,  =  a  component  equation  which  is  a  function  of  parameter  p4  only. 

C  =  constant 

Note  that  this  type  of  solution  is  often  sought  in  the  solution  of  linear  partial  differential 
equations  and  is  known  as  separation  of  variables.  For  example,  in  the  case  of  a  circular 
plate  placed  on  a  linear  elastic  medium  and  subjected  to  a  uniformly  distributed  load,  the 

real  response,  is  given  on  page  350  of  Das[47]  as 


0=  1.58 


where 


(5) 


Q  =  the  plate  sinkage. 

E  =  Young's  Modulus,  E=pj. 

\)  =  Poisson's  ratio,  \>=p2. 

q,b  =  constants  (respectively,  uniformly  distributed  load  and  plate  radius). 

Equation  (5)  is  a  multiplication  of  two  functions  of  the  parameters,  pj=E  and  P2=t>,  i.e., 

fj=l/E  and  f2=(l-\>2).  Therefore,  in  this  case  the  response  surface,  F  can  be  represented 
by  a  multiplication  of  the  component  equations  as  we  assumed  in  equation  (4).  However, 
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in  general  such  a  representation  is  accurate  only  in  a  small  region  due  to  geometric  and/or 
materia]  nonlinearities  in  the  system.  The  error  is  expected  to  be  small  if  the  range  of  p*  is 
small  for  each  of  the  m  parameters. 

Thus  the  process  of  building  the  response  surface  requires  holding  all  relevant  factors 
except  parameter  pj  constant(i.e.  geometry,  loading,  all  other  material  properties  pjt  j=l,2 
....  m  but  j&)  and  determining  the  component  equation  f(pj).  Once  all  the  component 
equations  are  determined,  equation  (4)  can  be  used  to  build  the  response  surface.  It  should 
be  recognized  that  for  each  given  geometry  and  loading  there  will  be  one  response  surface. 
In  the  case  of  plate  sinkage  tests,  for  a  given  plate  size  and  load  level  there  will  be  a 
response  surface.  Since  there  are  m  unknown  parameters,  at  least  m  field  measurements 
are  needed  to  solve  for  these  m  parameters.  In  practice,  it  is  preferable  to  have  more  than  m 
points(  i.e.  n>m)  so  that  the  m  parameters  can  be  determined  with  the  help  of  an 
optimization  algorithm.  Since  each  unreplicated  in-situ  measurement  corresponds  to  a 
given  geometry  and  loading,  each  of  these  experimental  values  correspond  to  a  point(  or 
contour)  on  one  response  surface.  Thus  each  of  the  n  unreplicated  measurements  will 
correspond  to  a  point(  or  contour)  on  one  of  the  n  distinct  response  surfaces.  Note  that 
more  than  one  observations  at  a  given  geometry  and  loading  refer  to  the  same  point  (or 
contour)on  a  response  surface  that  corresponds  to  that  geometry  and  loading.  Thus 
replicates  do  not  provide  additional  equations  to  solve  for  the  parameters,  but  help  in 
controlling  experimental  error.  Upadhyaya  et  al.[48)  suggested  that  at  least  eight  replicates 
to  adequately  deal  with  the  spatial  variability  in  the  case  of  in-situ  plate  tests.  Suppose  we 
have  n  distinct  combination  of  geometry  and  load  level  there  will  be  n  response  surfaces, 
Fj,  i=  1,2 n.  From  equation  (4)  we  get 

Ft  =  C1f11f12--.fiIn 
F2  *  C2f21f22—f2m 

(6) 


~  Cnfn|fn2...fnin 

where  fy  is  the  component  equation  corresponding  to  response  i  and  parameter  pj  and  Q  is 
the  constant  corresponding  to  the  same  response  surface  i. 

Since  each  of  the  material  parameter  has  its  own  range,  some  properties  such  as  Poisson's 
ratio,  u  vary  in  a  very  narrow  range  (0.0  to  0.5)  whereas  others  such  as  Young’s 
modulus,  E  may  vary  over  a  very  large  range(thousands  of  kPa).  From  the  point  of 
optimization  as  well  as  orthogonal  regression,  it  is  preferable  to  map  each  of  the  parameter 
to  the  same  range  through  scaling[32,49].  Each  of  the  unknown  parameter  was 
nondimensionalized  and  mapped  to  vary  from  -1  to  +1  by  the  following  transformation: 


P 


2<Pi  -  Pi> 

Pi  max  —  Pi  min 


(7) 


where: 

Pi  a  nondimensional  value  of  parameter  i. 
Pi  =  mid  point  value  of  parameter  i. 

Pi  max  =  upper  bound  value  of  parameter  i. 

Pi  min  -  lower  bound  value  of  parameter  i. 
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The  value  of  the  mid  point  is  zero,  upper  bound  is  1  and  the  lower  bound  is  - 1  for  each  of 
the  nondimensionalized  parameter. 


Let  fy  be  the  nondimensionalized  component  equation  corresponding  to  the 
nondimensionalized  parameter  p'j  and  test  condition  i.  The  relationship  between  fy  and 
fy  is  given  by 


where 


(8) 


Fj  as  computed  value  of  the  response  surface  F-  for  test  i  when  all  the  parameters  are  set 
equal  to  the  mid  point  value  of  zero. 


Moreover,  it  is  convenient  if  we  nondimensionalize  the  system  response  to  avoid  numerical 
problems  in  the  analysis.  The  nondimensionalized  response  surface  is  given  by: 

F'i  =  C,ifilfi2...fim  i=  1.2,3 . n  (9) 

Where 

F'i  =  nondimensionalized  response  surface  values  corresponding  to  the  i*  test  condition, 
C'j  =  correction  constant,  approximately  equal  to  1. 

The  data  for  the  creation  of  response  surfaces  can  be  obtained  from  any  analytical  or 
numerical  models.  We  propose  to  use  an  orthogonal  regression  technique  to  determine  the 
component  equation  fy.  The  use  of  an  orthogonal  regression  technique  not  only  provides 
an  equation  to  accurately  predict  the  overall  system  response,  but  also  provides  an  accurate 
estimate  of  regression  parameters[49,50,51].  An  accurate  estimation  of  regression 
parameters  is  essential  in  order  to  identify  the  unknown  material  parameters  by 
optimization.  The  function  f  y  is  an  orthogonal  polynomial  of  parameter  p  j  and  is  given 
by:  J 

k 

fij  =  X  (10) 

r=0 

The  values  of  ay.  ,r=l,2 . k  are  determined  by  using  model  response  (analytical  or 

numerical  such  as  FEM)  and  orthogonal  regression  techniques.  Only  requirement  for  the 
use  of  orthogonal  regression  in  curve  fitting  is  that  p  j  be  equally  spaced  during  model 
evaluation  while  all  other  material  parameters  be  held  at  the  mid  point  values.  The 
theoretical  value  of  the  correction  constant,  C'i  in  equation(9)  is  one.  However,  when 
curve  fitting  is  employed  to  determine  the  regression  coefficients^,  the  value  of  this 
correction  constant  may  be  slightly  different  than  one.  The  actual  value  of  C\  can  be  found 
by  a  employing  linear  regression  technique  between  Fj  and  (fiif  a— f  im)-  To  accomplish 
this  linear  regression,  model  response  at  orthogonal  points  used  in  building  the  response 
surface  and  some  additional  random  points  may  be  used. 

Higher  Order  Correction: 


As  stated  previously,  in  general  the  orthogonal  response  surface  is  expected  to  be  close  to 
the  true  model  response  only  near  the  mid  point  and  the  parameter  axes(p'i  axis).  As  we 
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start  moving  away  from  the  origin  or  the  parameter  axes,  the  two  surfaces  will  depart  from 
each  other.  At  large  distances  from  the  origin  and  the  parameter  axes  this  error  can  be 

significant  The  relation  between  the  nondimensionalized  true  response,  O'j  and  the 
response  surface,  Fj  is  given  by: 

+  (11) 
where  e,  is  the  error  in  our  approximation  ,<D'j  =  F,. 


By  assuming  that  the  function  <£ ,  is  "well  behaved"(  i.e.  analytic  everywhere  in  the 
parameter  space), this  function  can  be  represented  in  a  Taylor  series  as  follows: 


<D 


=  1  +  b,p',  +  b2p'2  +  ...  +  bmp'm  +  bjjpv  +  b^p'jP'2  +...+  bimP'lP’m  +binp’|3 
+  bU2p',2P'2  ♦•••+  bun^iy*  +  bjjjP’iP'jP’a  +...+  bnmP'lP'zP’ffl  +•••  02) 


where  coefficients, bj,  by,  by*.  etc.  for  i=l,2,  ...m;  j=l,2  ...,m;  k=l,2, . m  are 

respectively  related  to  the  partial  derivatives  of  the  function,  O'j  with  respect  to  p*i,  p'ip  j, 
P'i  Pj  P'k  etc.  at  the  origin(mid  point).  Equation  (12)  reduces  to  fy  along  p  j  axis.  i.e. 


bjPj  +  bjjPj2  +  bjjjpj3 


(13) 


Using  equations  (9),  (1 1),  (12)  and  (13)  we  get. 


=  dl2P*lP’2  +•••+  ^tmPlP'm  +  ^P'lP'l  +•••+  +••• 

+dl12P,l2P2  +-+  <JllmP,12P’m+  ^P'lPlP's  +- 


(14) 


where  "d"  s  are  constant  coefficients  related  to  the  cross  derivatives  of  <b'i  at  the  origin.  It 
should  be  noted  that  strictly  from  a  theoretical  point  of  view,  an  orthogonal  response 
surface  can  be  created  based  on  equation  (12)  rather  than  equation  (9)  which  relies  on  the 
product  of  component  equations.  In  such  a  case,  very  little  difference  is  expected  between 
the  real  surface  and  orthogonal  response  surface.  If  nine  equidistant  values  of  each  of  the 

parameter  p'i,  i=l,2 . on  are  used  in  evaluating  real  surface,  9™  model  evaluations 

will  be  needed.  If  m=2  then  81  model  evaluations  are  needed.  On  the  other  hand,  if  m=6, 
then  an  astronomical  531441  model  evaluations  are  necessary.  In  most  real  problems, 
where  FEM  evaluation  of  a  complex  model  is  necessary,  using  equation  (12)  as  a  basis  for 
the  response  surface  is  infeasible  except  for  the  case  of  a  two  parameter  model.  The 
response  surface  represented  by  equation  (9)  requires  only  [8*m+l]  model  evaluations( 
i.e.  for  m=2, 17  model  evaluations  are  necessary  whereas  for  m=6, 49  model  evaluations 
are  necessary). 


Second  Order  Correction: 


In  practice,  equation  (14)  will  be  truncated  at  some  convenient  point  The  truncated 
function  is  an  approximation  to  Ci  and  is  called  the  correction  function.Ej .  If  we  limit 

ourselves  to  only  the  product  of  the  type  p’jp 'j  for  i=l,2,  . m  and  j=l,2,  ....,m,  but 

i^j,  then  E,  will  be  a  second  order  function.  This  second  order  function,  Ej  contains  nc  min 
unknowns  given  by : 
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m(m-l) 

"c”"=  2  (15) 

In  order  to  determine  the  second  order  correction  function,  Ei  model  responses  are  obtained 
at  nc  additional  check  points,  where  tic  is  greater  or  equal  to  nc  min.  The  additional  check 
points  can  be  selected  randomly  or  in  a  deterministic  way.  It  can  be  shown  that  the  form  of 
the  second  order  correction  is: 


m-l  m 

-XX 

j=l  k=j+l 


e  (j-l)(2m-j)  .  P  jP'k 


k-j 


(16) 


The  "e"  coefficients  can  be  derived  from  a  set  of  nc  linear  equations  with  nc  ^  unknowns. 
A  multiple  linear  regression  technique  based  on  equation(16)  can  be  used  to  estimate  the 
"e"  coefficients.  Modification  to  the  response  surfaces  can  be  accomplished  by  adding 
equation  (16)  to  equation  (9).  The  resulting  improved  response  surface  is  given  by: 

Fi  =  Qfiifi2...fin,  +  Ei  i=l,2,3,...,n  (17) 

It  is  important  to  emphasize  that  the  second  order  correction  neglects  all  higher  orders  of 

£i-  There  may  be  some  situations  where  these  higher  order  corrections  are  necessary.  In 
such  cases,  it  is  possible  that  the  second  order  correction  may  even  give  poorer  results  than 
not  including  any  corrections.  In  situations  like  these,  use  of  equation  (9)  may  give 
more  accurate  results  than  equation  (17).  More  discussion  on  this  important  issue  will 
follow  when  we  consider  examples. 


Third  Order  Correction: 


In  order  to  get  more  accurate  results  to  the  function  E;  we  should  consider  the  higher 
order  corrections.  In  this  study  we  will  limit  ourselves  to  a  third  order  correction.  The  third 
order  correction  consists  of  all  cross  product  terms  of  the  parameters  upto  and  including 
the  third  order  terms.  The  third  order  function  Ej  is  the  summation  of  n^in  combination  of 
cross  products,  therefore  we  have  nc  ^  unknown  coefficients.  It  can  be  shown  that  the 
number  of  combinations,  nc  ^  is  given  by: 


_  3m(m  - 1)  1 

ncmin  —  -  +T 


(m-l)(m-2)  m-2  2 

<5  +  2-  J 

2  j=l 


The  function  Ej  for  the  third  order  correction  is: 


(18) 


Ej=mi  I  eiia(mJik)PjP’k  +  ni2  “l  I  eip(m  jkl)PjPkpi  + 

j=l  k=j+l  j=l  k=j+l  l=k+l 

2  £ei.,(n.j.k)Pj  Pkd-5)k) 

j=l  k=l 


(19) 


Where: 
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j  =  k 

j*  k 


5*=C 


2 

P(ra.  j.k.1)  =  +  i|i^t2m(m  -  j  - 1) + j] + Sr2  J  ■ 


iLlill 


(2m-k+j-2)  +  l-k 


(20) 

(21) 


(22) 


Y(m.j.k) 


m(m-l)  f 
2 


1  |~(m-l)(m-2) 

2|_  2 


+  °I2r2  +  (m-l)(j-l)  + 

r=l  . 


6 


(23) 


5  = 


k<j 

k>j 


(24) 


Once  again,  "e"  coefficients  can  be  derived  from  a  set  of  nc  (nc£nc  linear  equations 
with  nc  aii0  unknowns.  A  multiple  linear  regression  technique  based  on  equation  (19)  can 
be  used  to  estimate  the  "e"  coefficients.  The  modified  response  surface  is  given  by  equation 


(17). 


Estimation  of  Material  Parameters: 

Let  U4  be  one  of  the  n  independent  experimental  observations.  In  order  to  make  Uj 
consistent  with  F’i,  we  transform  it  into  a  nondimensional  value,  U'j.  The  relation 
between  Ui  and  U'i  is  given  by: 

U’j  =  ^  (25) 

The  subtraction  of  equation  (25)  from  equation  (17)  yields  a  set  of  n  nonlinear  equations  in 
m  unknowns. 


Fi  *  U*!  =  0 
F2  -  U’2  =  0 

(26) 


F„  -  U’n  =  0 

In  general  equation  (26)  is  seldom  an  equality  due  to  the  presence  of  approximation  as  well 
as  experimental  errors.  One  method  of  determining  engineering  properties  of  the  material 
involves  optimizing  sum  of  squares  of  residuals,  SSR  defined  by: 

>  ' 

Z<Fi  ~  U'.)2  ’ 

.  pI 


SSR  =  min  < 


(27) 
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The  expression  for  the  SSR  in  equation(27)  is  used  as  an  objective  function  and  a  nonlinear 
optimization  technique  is  used  to  solve  for  material  parameters.  Since  F'i  is  an 

approximation  to  and  Fj  can  be  obtained  from  equation  (9)  which  corresponds  to  the 
original  orthogonal  response  surface  or  equation  (17)  which  includes  higher  order 
correction,  we  can  get  different  solutions  in  the  vicinity  of  the  real  solution  depending  on 
which  expression  for  Fj  is  used.  The  SSR  of  these  solutions  is  of  the  same  order,  thus 
the  minimum  value  of  the  SSR  does  not  necessarily  indicate  the  best  solution.  Of  course, 

the  best  solution  can  be  obtained  from  the  sum  squares  of  the  residuals  which  uses  d>'j  in 
equation  (27),  i.e. 


SSR*  =  min  i 


!>',  -  u\r 


l  i=l 


(28) 


Where  SSR*  is  the  minimum  residuals  of  the  real  response  surfaces  and  the  test  results.  It 
is  recommended  to  use  equation  (28)  at  these  different  optimum  solutions  suggested  by 
different  versions  of  Fj  to  obtain  the  best  results.  We  will  return  to  this  question  when 
we  consider  example  problems. 

Transformation  of  Experimental  Results: 

As  explained  previously,  the  accuracy  of  function  Fj  may  be  high  in  some  region  on  the 
response  surface  and  low  in  another  region.  This  implies  that  the  estimated  parameters 
will  be  high  in  accuracy  sometimes  and  poor  in  accuracy  some  other  times.  Generally 
speaking,  the  inaccuracy  increases  as  we  move  away  from  the  origin  and  parameter  axes. 
At  the  origin,  the  nondimensionalized  F’i  has  a  value  of  unity[cf.  eq.  (9)].  In  some  sense, 
as  the  values  of  F'i  change  from  unity  the  difference  between  the  real  and  the  response 
surface  values(  both  corrected  and  uncorrected)  tend  to  increase.  Thus  the  values  of  Fj 
can  be  used  as  a  measure  of  this  departure.  This  argument  suggests  that  there  exists  a 

function  <J>'i  =<J>’i(F’i).  Inverse  of  this  transformation,  F|=  Fj(<I>'i)  is  of  particular 
interest  in  our  case.  This  relationship  can  be  used  as  a  transformation  rule  for  experimental 

values  by  replacing  the  real  response,  O'i  by  experimental  value,  U’i.  If  we  denote  the 
transformed  experimental  value  which  corresponds  to  F'i  by  U*j,  then  we  have 
U*i=U*j(U'i).  The  transformation  function  can  be  obtained  by  conducting  a  polynomial 

regression  between  Fi  and  <3>'i  in  some  acceptable  range,  Oj  min  and  <J>j  max ,  thus: 


(29) 


where  "gj"  s  are  regression  coefficients.  The  corresponding  transformation  rule  for  the 
experimental  values  is  given  by: 


u‘  =  X  SjU'i 

j=0 


(30) 
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The  value  of  U*j  calculated  from  equation  (30)  can  be  used  for  replacing  U*  in  equation 
(26)  and  (27).  This  modification  may  significantly  improve  the  accuracy  of  the  parameters 
obtained  through  optimization.  We  will  explore  this  aspect  in  more  detail  when  we  consider 
the  example  problems. 

COMPUTER  IMPLEMENTATION  OF  THE  MODEL: 

A  FORTRAN  program  was  developed  to  implement  this  inverse  solution  technique  to 
estimate  material  parameter  values.  This  task  involves  several  steps  as  shown  in  figure  1. 
An  interface  program  is  used  to  transform  the  analytical  or  numerical  results  to  a  format 
accepatable  to  our  inverse  solution  program.  The  inverse  solution  program  uses  model 
response  (either  analytical  or  numerical  such  as  FEM  analysis)  for  the  creation  of  the 
orthogonal  response  surface.  Higher  order  correction  are  implemented  on  to  this 
orthogonal  response  surface  in  the  next  step.  Following  this  expeimental  data  are  input 
and  a  transformation  is  performed  on  these  data  to  relate  them  to  response  surface  points. 
In  the  last  step,  an  optimization  procedure  is  carried  out  to  obtain  the  best  estimates  of  the 
values  of  material  parameters.  This  program  allows  for  selection  of  any  one  of  the 
following  optimization  subroutines  called  from  IMSL  library: 

1.  Nonlinear  least  squares  techniques, 

2.  Complex  algorithm, 

3.  Quasi-Newton  method, 

4.  Modified  Newton  Method, 

5.  Conjugate  Gradient  Method. 

CASE  STUDIES: 

Two  Parameter  Model: 

A  simple  two  parameter  material  model  is  selected  to  illustrate  the  main  features  of  this 
technique.  A  cylindrical  bar  or  soil  column  under  uniaxial  compression  is  considered.  The 
bar  is  made  of  a  hypoelastic  material  with  an  incremental  constitutive  law  given  by(p 
139,Desai  and  Siriwardane  [17]): 

do  =  (Eq  +  EjCride  (31) 

where  Eq  and  Et  are  the  material  parameters  of  interest  in  this  model,  o  is  axial 
stress(positive  in  compression),  and  e  is  axial  strain.  Integration  of  equation  (31)  yields: 


where: 

u  =  deformation  of  the  bar  (the  contraction). 

L  =  length  of  the  undeformed  bar. 

Equation  (32)  will  be  used  as  a  basis  to  build  the  response  surface,  Fi.  The  response 
surface  will  be  developed  in  the  following  range  of  parameters  Eq  and  Ej: 

E0  ^  689.5  kPa  (100  psi)  Eq  max=  6894.8  kPa  (1000  psi) 

Elnun=  10.0  Elmax=  100.0 

E0oudPoin«=3792.1  kPa  (550  psi)  E! 

mid  point2 


55.0 


13 


The  nondimensional  parameters  (equation  7)  for  this  case  are: 


2(Eq  Eq  jjjjd  point) 


2(Ej  Ej  nuj  point) 


Eo  =  ~p - _  p  and  E,  -  - - 

R)  max  *-Omin  C1  mu  **1  mil 

The  response  at  the  mid  point  of  the  parameters  is(i.e.  origin): 


— taf 

*-’1  mid  point  V 


EQ  mid  point  ^1  mid  point*- 
Eq  mid  point 


(33) 


(34) 


The  nondimensional  representation  of  the  real  surface  is  obtained  by  dividing  equation  (32) 
by  equation  (34).  The  plot  of  the  nondimensionalized  real  surface  obtained  by  using  an 

applied  stress  of  a=  689.5  kPa(100  psi),  is  given  in  Fig.  2.  The  approximation  of  the 
surface  without  any  higher  order  correction  is  shown  in  Fig.  3,  and  the  error  of  this 
approximation  is  shown  in  Fig.  4.  The  response  surface  describes  the  real  function  with 
reasonable  accuracy  except  at  the  comers.  The  error  is  particularly  high  as  both 
parameters  (Eq  and  E,)  approach  their  minimum  values.  The  response  surface  for  this  case 
which  includes  second  order  correction  is: 

Fi  =  f^E’olfjfE1!)  +  ei  E'qE’j  (35) 

where  f !  and  f2  are  orthogonal  polynomial  functions  of  E'0  and  E',  respectively  and  e4 
is  second  order  correction  coefficient .  The  function  f  j,  f2  and  the  coefficient  ct  were 
found  as  described  previously. 

The  second  order  correction  for  this  case  was  obtained  by  using  the  edge  points,  the  mid 
points  of  the  lower  and  upper  range  for  each  parameter  -  a  total  of  16  combinations.  The 
response  surface  with  second  order  correction  is  shown  in  Fig.  5,  and  the  difference 
between  the  real  values  and  the  response  surface  values  is  plotted  in  Fig.  6.  The  second 
order  correction  decreases  the  error  in  the  zones  of  high  error  (e.g.  in  the  region  near  the 
minimum  values  of  Eq  and  Ei).  However,  this  correction  to  the  response  surface  increased 
the  error  in  some  other  regions  where  the  error  was  negligible  previously.  The  basic 
assumption  of  the  second  order  correction  is  that  all  higher  order(third  order  and  higher) 
cross  products  are  negligible.  In  this  particular  example,  the  plot  of  the  error  shown  in  Fig. 
4  indicates  a  high  curvature  in  a  small  region  [in  fact  only  about  3%  region]  where  both 
parameters  approach  their  minimum  values.  Thus,  second  order  correction  does  not 
represent  the  error  property  everywhere  in  the  region. 

Effect  of  including  the  third  order  correction  (TOC)  was  also  examined  for  thi.»  two 
parameter  case.  The  form  of  the  third  order  correction  is: 

2  2  (36) 

Ej  —  ejjE()Ej  +  ei2E0Ej  +  C13E0E1 

Figure  7  shows  the  effect  of  including  third  order  correction  on  the  response  surface. 
Inclusion  of  third  order  correction  further  reduces  the  error  in  the  zone  of  high  error  (i.e.  in 
the  3%  region  where  Eq  and  Ei  values  are  near  or  at  their  minimum).  Figure  8  is  a  plot  of 


14 


error  when  the  third  order  correction  is  included.  Comparison  of  this  figure  with  figures  4 
and  6  shows  that  although  the  response  surface  which  includes  third  order  correction 
reduces  the  error  in  the  zone  of  high  error ,  the  error  in  other  regions  does  not  necessarily 
decrease.  In  fact,  the  error  increases  slightly  in  some  areas. 


In  order  to  determine  the  material  parameters(Eo  and  Ei),  seven  different  resoponse 
surfaces  were  created  using  seven  different  applied  stresses.  The  stress  values  used  were 
344.7  kPa  (50  psi),  517.1  kPa  (75  psi),  689.5  kPa  (100  psi),  861.8  kPa  (125  psi), 
1034.2  kPa  (150  psi).  11206.6  kPa  (175  psi)  anu  1379.0  kPa  (200  psi).  The  adequacy 
of  the  method  is  illustrated  by  12  examples.  Table  1  lists  the  parameter  values  selected  for 
the  simulation  purpose.  The  first  set  of  parameter  values  are  randomly  selected.  The  other 
eleven  examples  are  based  on  parameter  values  along  the  diagonal,  E'o=E'i.  The  equation 
(32)  was  used  to  calculate  the  real  response.  These  values  were  used  instead  of  the 
experimental  values  in  the  optimization  step.  Since  there  is  no  experimental  error  in  this 
case,  we  should,  in  principle,  get  exact  values  of  the  assumed  parameters  back.  Inaccuracy 
in  the  results  is  solely  due  to  die  inadequacy  of  the  response  surface.  Five  different  initial 
guess  values  of  the  parameters  were  considered  in  the  nonlinear  optimization  process  for 
each  one  of  these  examples.  The  first  initial  guess  values  were  the  mid  point  values  of  the 
parameters,  three  others  were  selected  randomly  and  the  fifth  one  was  the  exact  solution. 

The  transformation  equation  was  estimated  using  50  random  points  of  the  real  surface 
using  equation  (29)  as  a  basis  fv,  regression. 

Example  1:  Eq  and  Ej  were  selected  randomly: 

The  values  of  the  simulated  parameters  were  :  Eq=  5666.1  kPa  (821.8  psi)  and  E!=61.8. 
The  nondimensionalized  real  function  value  at  an  applied  pressure  of  689.5  kPa  (100  psi) 
is  0.80.  The  results  of  optimization  are  listed  in  Table  2. 

Note  that  both  uncorrected  response  surface  and  respose  surface  with  third  order  correction 
resulted  very  good  solutions  with  negligible  errors.  However,  the  second  order  correction 
yielded  relatively  poorer  results.  When  the  transformation  based  on  equation(30)  was 
employed,  all  three  correction  methods  yielded  reasonably  good  results  (Table  3). 
Transformation  technique  significantly  improved  the  parameter  estimation  for  the  case  in 
which  no  correction  was  employed.  This  transformation  was  beneficial  in  reducing  error 
for  the  second  order  correction  technique  also,  particularly  at  low  parameter  values. 
However,  the  transformation  technique  was  not  quite  as  beneficial  in  the  case  of  the  third 
order  correction  method.  Even  in  this  case  the  error  in  estimation  of  the  parameters  were 
reduced  slightly.  Thus,  in  general  transformation  technique  leads  to  more  accurate 
parameter  estimation. 

It  should  be  noted  that  this  particular  soil  model  shows  large  increases  in  response  when 
Eo  and  Ei  values  are  very  small  compared  to  all  other  values  of  Eo  and  Ej.  This  portion 
of  the  graph  corresponds  to  only  3%  of  the  parameter  range  (low  values  of  the 
parameters).  This  makes  it  difficult  to  generate  a  response  surface  which  is  good 
everywhere  in  the  region.  Higher  order  corrections  tend  to  predict  the  response  better  in 
this  region.  In  so  doing,  they  become  less  accurate  in  other  regions  -  especially  at  large 
values  of  the  parameters  (see  Figs,  3  through  8).  This  is  an  unusal  situation  which 
resulted  in  the  uncorrected  surface  to  generally  estimate  parameter  values  more  accurately 
than  when  the  third  order  correction  was  employed  along  with  the  transformation 
technique.  In  a  well  behaved  system  (i.e.  no  singularities  or  large  increases  in  responses 
for  small  changes  in  parameter  values)  the  method  which  employs  third  order  correction  is 
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expected  to  yield  more  accurate  parameter  valves.  In  fact,  even  very  complicated  models 
do  not  show  such  singularities  or  large  changes  in  a  small  region  as  we  will  see  with  a  five 
parameter  Drucker-Prager  model  described  below. 

A  Five  Parameter  Nonlinear  Elastic  Soil  Model  with  Extended  Druker 
Prager  Yield  Criteria: 

The  elasto-plastic  constitutive  material  model  with  Druker-Prager  yield  criteria  is  widely 
used  in  geomechanics[45,46].  We  assume  that  in-situ  tests  consist  of  plate  penetration 
tests  using  circular  plates.  Here  we  will  not  consider  the  actual  field  data.  Analysis  of  the 
field  data  to  identify  material  parameters  will  be  dealt  later.  We  explore  the  feasibility  of  the 
proposed  methodology  in  this  example  for  this  fairly  complex  material  model.  A 
commercial  finite  element  program,  ABAQUS  was  used  in  this  study  to  obtain  model 

response.  The  orthogonal  response  surfaces  were  built  using  six  parameters  k/d,K,C,<}> 
and  e. 

Two  plates  of  diameters  50  mm(2  in.)  and  100  mm  (4  in.)  were  simulated,  with  applied 
load  ranging  from  137.9  kPa(  20  psi)  to  1034.2  kPa  (150  psi)  in  increments  of  68.9  kPa( 
10  psi),  a  total  of  14  tests  for  each  plate.  A  response  surface  was  built  for  each  of  those 
tests  in  the  following  parameter  range: 
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A  typical  plot  of  the  non-dimensional  sinkage  as  a  function  of  non-dimensional  values  of 

parameter  k  for  a  100  mm  (4  in)  plate  subjected  to  551.6  kPa  pressure  (80  psi)  is  shown  in 
Fig.  9.  Note  that  all  other  parameters  are  held  constant  at  their  corresponding  midpoint 
values.  Figures  10  through  14  are  similar  plots  except  that  the  dependent  variable  has  been 

changed  to  o,  K,  C,  $,  and  e  respectively.  These  curves  show  an  extremely  good  fit 
between  the  real  response  curve  and  the  response  curve  obtained  by  the  orthogonal 
regression  (R2  >0.997  for  all  cases). 

The  graph  of  the  real  surface  versus  the  orthogonal  response  surface  without  any 
correction  for  a  100  mm  (4  in.)  plate  subjected  to  an  applied  load  of  103.4  kPa  (15  psi) 
is  shown  in  Fig.  15.  This  graph  consists  of  55  orthogonal  points  and  an  additional  60 
random  points.  Figure  16  is  similar  to  Fig.  15  except  that  the  applied  load  is  551.6  kPa 
(80  psi)  in  this  case.  When  the  applied  pressure  low,  soil  deformation  is  small  and  the 
soil  medium  behaves  similar  to  elastic  materia  m  fact,  in  our  case  soil  is  modeled  as 
nonlinear  elastic  but  if  displacements  are  small  even  nonlinear  behavior  can  be 
approximated  by  linear  behavior).  In  this  case  the  real  surface  and  orthogonal  surface  are 
almost  identical  (Fig.  15).  As  the  load  is  increased,  soil  will  yield  and  subsequent  plastic 
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flow  will  take  place  as  per  the  assumed  model.  Figure  16  reveals  that  under  high  load  the 
orthogonal  response  surface  begins  to  depart  from  the  real  surface  when  non- 
dimensionalized  displacements  are  below  0.5  or  exceed  1.5.  Figures  17  and  18  are  similar 
to  figures  15  and  16  except  that  a  second  order  correction  has  been  added  to  the  orthogonal 
response  surface.  These  figures  indicate  that  there  has  been  only  marginal  improvements  in 
these  curves  (especially  Fig.  18).  Perhaps  a  higher  order  correction  is  beneficial  especially 
at  high  plate  loads.  Figures  19  and  20  are  similar  to  figures  15  and  16  (also  17  and  18) 
respectively,  except  that  a  third  order  correction  has  also  been  added  to  the  orthogonal 
response  surface.  Inclusion  of  third  order  correction  has  resulted  in  an  orthogonal 
response  surface  which  is  almost  identical  to  the  real  surface  even  at  high  loads.  This 
indicates  that  an  orthogonal  response  surface  with  third  order  correction  can  be  used 
reliably  to  predict  the  real  response  without  having  to  resort  to  FEM  analysis. 


Parameter  Estimation: 


Since  we  are  dealing  with  a  nonlinear  problem  the  solution  is  not  necessarily  unique. 
Following  recommendations  may  be  used  as  a  guide  for  selecting  the  best  solution  from 
several  optimum  solutions  resulting  from  the  presence  of  "local  minim ums": 

1.  Discard  all  solutions  that  have  a  significantly  high  SSR  (cf.  equation  27). 

2.  Use  more  than  one  geometry  (i.e.  50  mm  and  100  mm  diameter  plates)  and  look  for 
optimum  for  each  of  the  geometries  and  also  the  combination  of  all  the  geometries.  Accept 
those  solutions  which  are  approximately  same  in  all  cases.  From  a  practical  point  of  view 
two  plates  will  be  sufficient 

3.  Reject  any  solution  in  which  more  than  one  parameters  hit  the  bounds  of  the  search 
domain.  The  probability  of  more  than  one  parameter  hitting  the  bounds  simultaneously  is 
low.  If  in  fact,  if  this  really  is  the  case  for  several  initial  guesses,  and  the  above  two  criteria 
will  be  met. 

4.  In  spite  of  these  steps,  if  more  than  one  optimal  solutions  are  obtained,  we  recommend 
the  use  of  equation  (27)  to  compute  SSR.  If  the  response  surface  with  higher  order 
corrections  has  been  properly  verified  as  good  (i.e.  using  figures  such  as  19  and  20 
described  earlier),  then  the  solution  which  yields  the  minimum  SSR  should  be  accepted  as 
the  best  solution.  If  possible  one  could  use  equation  (28)  to  compute  SSR*  to  provide 
additional  verification.  Although,  such  an  approach  is  preferable,  evaluation  of  equation 
(28)  requires  some  limited  FEM  analysis  (i.e.  one  for  each  plate  for  each  competing 
solution). 

To  explore  the  suitability  of  this  method  to  identify  the  material  parameters  of  this  complex 
constitutive  equation  for  soil,  we  selected  five  different  random  sets  of  parameters  and 
conducted  simulation  studies  to  obtain  true  response.  Subsequently,  these  true  responses 
were  used  as  inputs  into  to  the  response  surface  methodology  to  re-predict  those 
parameters.  Since  third  order  correction  (TOC)  to  the  orthogonal  response  surface  appears 
to  be  necessary  to  obtain  reasonable  results,  we  will  only  explore  the  situation  in  which 
TOC  is  added  to  the  orthogonal  response  surface.  We  chose  five  different  random  sets  of 
parameter  values  to  be  re-predicted  using  the  optimization  technique.  These  random  sets 
of  points  are  as  follows: 
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The  initial  guess  values  selected  were  the  "exact  solution",  "mid  point  values"  and  a  set  of 
five  randomly  selected  parameter  values  listed  below: 
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The  re-prediction  process  was  carried  out  using  50  mm  (2  in. )  plate,  100  mm  (4  in.)  plate 
and  a  combination  of  50  mm  (2  in.)  and  100  mm  (4  in.)  plates.  For  each  case  we  used  the 
exact  solution  and  the  mid  point  as  initial  guesses.  The  five  randomly  selected  guess  points 
were  used  only  for  point  #  1  for  both  plates  and  also  the  combination  of  plates.  However, 
for  the  100  mm  (4  in. )  plate  all  seven  initial  guess  values  were  use  to  seek  the  optimum 
solution  for  each  of  the  five  random  set  of  parameter  values. 

The  results  of  this  analysis  are  listed  in  Table  4.  An  examination  of  the  results  indicates 
that  the  reasonable  solution  with  minimum  SSE  usually  leads  to  good  solution  except  for 

point  #4.  Point  #4  results  in  very  large  errors  for  both  parameters  k  and  o.  However,  an 
examination  of  SSE  indicates  that  none  of  the  solutions  is  reasonable.  Our  suspicion  is 

that  for  this  values  of  Kando,  the  soil  is  extremely  hard  and  deforms  very  little.  Under 
these  circumstances,  the  nonlinear  elastic  model  for  soil  with  Drucker-Prager  yield  criteria 
is  perhaps  inappropriate. 

DETERMINATION  OF  IN-SITU  SOIL  PROPERTIES: 


Field  tests  were  conducted  during  November,  1991  and  September,  1992  using  our 
instrumented  soil-test  device  in  a  Yolo  loam  soil  in  the  vicinity  of  the  U.C.  Davis  campus. 
In  November  1991  tests  were  conducted  in  an  undisturbed  soil  using  50.8  mm  (2  in.), 
76.2  mm  (3  in.),  101.6  mm  (4  in.),  127  mm  (5  in.),  and  152.4  mm  (6  in.)  sinkage  plates. 
When  we  conducted  the  field  tests  we  were  under  the  impression  that  more  than  one 
geometry  (plate  sizes)  were  necessary  to  obtain  the  engineering  properties  of  soil  by 
theinverse  solution  technique.  However,  later  we  found  that  this  is  not  necessarily  the 
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case.  Even  use  of  one  plate  size  appears  to  be  sufficient  In  this  study,  field  test  results 
for  50.8  mm  (2")  and  101.6  mm  (4")  plates  were  only  used  to  determine  the  engineering 
properties  of  soil  by  the  inverse  solution  technique. 

Eight  replicates  were  obtained  for  each  plate.  Sinkage  test  data  were  analyzed  using 
Recee's  approach,  i.e. 

p  =  k  (z/r)®  (37) 


where 

p  -  applied  pressure 
k  =  sinkage  constant 
z  =  soil  sinkage 
r  =  plate  radius 
n  =  empirical  constant 

Table  9  lists  the  mean  values  of  sinkage  coefficients  for  each  of  the  plate  tested.  The 
values  of  k  and  n  for  the  50.8  mm  (2")  and  101.6  mm(4")  plate  were  used  in  estimating 
mean  field  response  corresponding  to  a  desired  pressure  for  a  given  plate  during  the 
optimization  process  to  "identify "  soil  parameters. 

Soil  shear  tests  were  conducted  using  two  different  grouser  plates  [Plate  #1:  203  mm  long 
x  76  mm  wide,  and  Plate  #2:  178  mm  long  x  86  mm  wide].  Each  plate  was  tested  at  two 
different  vertical  loads  and  each  test  was  replicated  three  times.  Grouser  plate  test  results 
were  analyzed  using  the  following  equation: 


j 

x  =  [c + p * tan(<|>)](l - e  K)  (38) 


where 

x  =  shear  stress,  kPa 

c  =  cohesion,  kPa 

p  =  pressure  on  the  plate,  kPa 

<j>  -  soil  internal  friction  angle 
j  =  shear  deformation,  mm 
K  =  shear  modulus,  mm 

A  nonlinear  regression  technique  was  employed  to  fit  the  data  to  equation  (38)  and  obtain 

shear  parameters.  Maximum  shear  stress,  Xma*  for  each  test  was  calculated  using  the 
following  equation: 

rmax=c  +  p*tan((j>)  (39) 

The  analysis  of  the  experimental  data  resulted  in  a  mean  value  of  cohesion  of  1 1.5  kPa 
and  soil  internal  friction  angle  of  32.9  deg.  for  these  tests. 

Cone  index,  bulk  density  and  moisture  content  data  were  also  obtained  in  the  test  site.  Eight 
replicates  of  cone  index  profiles  were  obtained  in  the  top  152.4  mm  (6")  layer.  The  cone 
index  values  were  averaged  over  the  depth  to  get  a  representative  cone  index  value  for  each 
location.  Subsequently,  all  eight  replicates  were  averaged  to  get  a  mean  cone  index  value 
for  this  particular  soil  condition.  Five  bulk  density  and  moisture  content  data  were  also 
obtained  in  the  test  site.  Average  cone  index  value  was  816  kPa,  dry  bulk  density  was 
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1510  kg/m3,  and  moisture  content  was  8.9%  (dry  basis).  The  void  ratio  was  0.755  based 
on  a  particle  density  of  2650  kg/m3. 

During  September,  1992  only  three  sinkage  plates  were  used  for  sinkage  tests.  Two 
distinct  soil  conditions  (  undisturbed  and  tilled/loose)  were  included  in  these  tests.  Once 
again  eight  replicates  of  sinkage  tests  were  obtained  for  each  plate  and  analyzed  using 
equation  (37).  Sinkage  parameters  for  the  undisturbed  and  tilled  soil  conditions  of  the 
November  1992  tests  are  also  listed  in  Table  5.  Once  again,  the  mean  sinkage  parameters 
corresponding  to  50.8  mm  (2")  plate  and  101.6  mm  (3")  plate  were  used  in  estimating  field 
response  for  identifying  Engineering  parameters  of  soil. 

Shear  test  procedure  as  well  as  data  analysis  techniques  were  similar  to  the  procedure 
employed  in  analyzing  the  November,  1991  shear  tests.  The  mean  value  of  cohesion  for 
the  undisturbed  soil  condition  was  32.3  kPa,  and  the  internal  angle  of  friction  was  27.2 
deg.  The  corresponding  values  for  the  tilled  soil  was  22.7  kPa  and  22.8  deg. 

Moreover,  soil  bulk  density,  cone  index  and  moisture  content  data  were  also  obtained  as 
described  for  the  November  1991  tests.  The  undisturbed  (also  referred  as  firm)  soil  had  a 
mean  dry  bulk  density  of  1510  kg/m3,  moisture  content  of  5.  14  %  (dry  basis),  and  an  void 
ratio  of  0.755.  The  loose  or  tilled  soil  had  a  dry  bulk  density  of  1433  kg/m3,  4.55  % 
moisture  content  (dry  basis),  and  a  void  ratio  of  0.851.  The  cone  index  data  were 
inconsistent  and  were  ignored  for  these  soil  conditions. 

Table  6  lists  the  Engineering  parameters  of  soil  estimated  from  the  optimization  process 
which  utilized  the  orthogonal  response  surface  including  the  third  order  correction.  Both 
the  best  results  based  on  SSE  and  SSR,  and  reasonable  results  based  on  our  search  criteria 
are  listed  in  Table  6  for  all  the  three  soil  conditions.  The  best  estimates  of  the  cohesion  and 
soil  internal  friction  angle  values  listed  in  Table  6  do  not  agree  with  the  corresponding 
values  listed  in  Table  5,  which  are  grouser  shear  test  results.  This  agrees  with  our 
hypothesis  that  the  grouser  shear  test  provides  geometry  dependent  soil  parameters,  but 
not  the  basic  soil  constitutive  property.  Only  the  101.6  mm  (4  in.)  plate  was  used  for 
"identifying"  soil  parameters  through  optimization.  Figures  21  and  22  show  the 
experimental  and  simulated  sinkage  for  a  101.6  mm  (4  in.)  plate  obtained  using  back- 
calculated  soil  parameters  for  September  1992  tests  in  an  undisturbed  soil  (firm  soil). 
These  results  indicate  that  the  estimated  soil  parameters  are  very  good.  However,  when 
these  same  parameters  were  used  to  compare  the  response  of  a  50.8  mm  (2  in.)  plate  in  the 
same  soil  condition  poor  agreement  was  found  between  experimental  and  simulated  sinkage 
( Figs.  23  ).  This  plot  indicates  that  the  parameters  predicted  from  101.6  mm  (4  in.)  plate 
tests  are  unable  (under  predict)  to  predict  the  behavior  of  50.8  mm  (2  in.)  plate  in  the  field. 
Similar  results  were  obtained  in  the  other  soil  conditions  tested  also.  We  feel  that  this  is  due 
to  the  edge  effect  which  is  not  included  in  our  model.  Use  of  an  interface  element  at  the 
soil-plate  interface  appears  to  be  necessary.  Since  101.6  mm  (4  in.)  plate  is  less 
susceptible  to  edge  effect  compared  to  the  50.8  mm  (2  in.)  plate,  we  feel  the  parameters 
estimated  from  a  101.6  mm  (4  in.)  plate  are  more  reliable.  We  recommend  using  a  large 
diameter  plates  in  futv  re  tests. 

CONCLUSIONS 

Based  on  this  study  we  reached  the  following  conclusions: 

1)  A  response  surface  methodology  based  on  an  orthogonal  regression  in  the  parameter 
space  has  been  developed  to  "identify"  ,  or  "calibrate"  engineering  properties  of  any 
material  based  on  in-situ  tests.  The  orthogonal  response  surface  was  created  from  an 
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analytical  or  numerical(such  as  FEM)  solution  to  the  underlying  differential  equation  of  the 
system  which  utilizes  these  engineering  properties  in  a  constitutive  equation.  A 
transformation  technique  was  developed  to  map  the  model  response  or  experimental  data  on 
to  the  response  surface. 

2)  The  proposed  methodology  worked  very  well  (i.e.  very  little  error)  in  the  case  of  a  two 
parameter  hypo-elastic  model  for  soil.  When  the  second  order  correction  was  included  with 
a  transformation  of  data  very  small  errors  resulted  in  parameter  estimation.  Inclusion  of 
third  order  correction  to  the  orthogonal  response  surface  reduced  the  chance  of  large  error 
in  parameter  values. 

3)  When  this  technique  was  used  in  the  presence  of  random  noise,  the  predicted 
parameters  were  found  to  be  insensitive  to  the  noise. 

4)  When  this  methodlogy  was  applied  to  a  complex  five  parameter  model  for  soil 
(nonlinear  elastic  behavior  with  Drucker-Prager  yield  criteria  and  associated  plastic  flow 
upon  yield),  it  appeared  to  work  reasonably  well.  A  third  order  correction  to  the 
orthogonal  response  surface  appears  to  be  necessary  to  obtain  reasonably  good  solution. 

When  both  the  logarithmic  bulk  modulus  (k)  and  Poisson's  ratio  (u)  are  low,  soil  becomes 
very  rigid  and  the  methodlogy  will  not  yield  a  good  solution.  Under  such  circumstances, 
perhaps  the  soil  model  chosen  is  inappropriate. 

5)  The  response  surface  methodology  was  successfully  employed  to  "identify" 
engineering  properties  of  soil  based  on  field  tests  for  different  soil  conditions  in  a  Yolo 
loam  soil.  We  suspect  that  edge  effect  makes  the  parameter  prediction  using  field  data 
corresponding  to  small  plates  such  as  50.8  mm  (2  in.)  diameter  plate  inaccurate.  Use  of 
larger  plates  such  as  101.6  mm  (4  in.)  plate  is  recommended  to  reduce  this  edge  effect 
Use  of  Teflon  coated  plate  with  beveled  edges  and  slip  elements  at  the  plate  edge  in  the 
model  may  also  increase  the  accuracy  of  parameter  prediction. 

The  proposed  methodology  has  not  only  applications  in  geomechanics,  but  also  in  other 
areas  such  as  biological  engineering(plants  and  animal  tissues,  food  products  etc.)  where 
non-destructive  in-situ  tests  are  the  only  means  of  obtaining  accurate  estimate  of 
engineering  parameters. 
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Table  1.  Parameter  values  used  in  the  simulation  studies 


Sequence 

number 

True  parameters  values 

Nondimensional  parameters 

values 

Eo(kPa) 

E, 

E'o 

E'i 

1 

5666.1 

61.8 

0.604 

0.151 

2 

999.7 

14.5 

-0.9 

-0.9 

3 

1310.0 

19.0 

-0.8 

-0.8 

4 

1930.5 

28.0 

-0.6 

-0.6 

5 

2551.0 

37.0 

-0.4 

-0.4 

6 

3171.6 

46.0 

-0.2 

-0.2 

7 

3792.1 

55.0 

0.0 

0.0 

8 

4412.6 

64.0 

0.2 

0.2 

9 

5033.2 

73.0 

0.4 

0.4 

10 

5653.7 

82.0 

0.6 

0.6 

11 

6274.2 

91.0 

0.8 

0.8 

12 

6894.8 

100.0 

1.0 

1.0 
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Table  2.  Parameter  values  obtained  through  optimization  for  the  cases  when  no  correction 
was  included  as  well  as  the  case  for  which  second  order  and  third  order  correction  were 

included 


%  % 
error  SSE  error 
ofE,  xlO3  ofEo 


1  Third  order  correction 

% 

error 
of  Eq 

% 

error 
of  Ej 

SSE 

XlO-3 

0.95 

1.48 

0.3104 

1.88 

4.32 

2.2190 

17.13 

3.30 

30.417 

17.02 

2.21 

29.463 

8.92 

1.12 

8.0740 

2.38 

0.37 

0.5784 

0.03 

0.00 

0.0001 

1.85 

0.44 

0.3632 

1.01 

1.74 

0.4030 

3.06 

0.15 

0.9370 

15.69 

8.03 

31.064 

25.21 

0.08 

63.575 
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Table  3.  Parameter  values  obtained  through  optimization  when  transformation  technique 

was  employed 


Without  correction 

Second 

order  correction 

Third  order  correction  \ 

No. 

% 

error 
of  Hq 

% 

error 
of  Ej 

SSE 
x  10-3 

% 

error 
of  Eq 

% 

error 
of  E, 

SSE 
x  10-3 

% 

error 
Of  Eq 

% 

error 
of  Et 

SSE 

X  10-3 

1 

0.86 

0.24 

0.0798 

1.98 

0.33 

0.4030 

1.78 

0.55 

0.3456 

4.09 

2.83 

2.4780 

2.81 

1.88 

1.1470 

0.56 

0.09 

0.0395 

B 

1.49 

3.92 

1.7567 

2.79 

0.09 

0.7796 

2.06 

5.47 

3.4132 

5.48 

2.83 

3.8048 

0.06 

1.75 

0.3060 

9.12 

4.33 

10.202 

5 

4.59 

2.39 

2.6734 

2.47 

2.62 

1.3003 

5.98 

2.49 

4.1949 

6 

4.14 

1.49 

1.9347 

2.01 

0.71 

0.4562 

3.95 

0.78 

1.6195 

D 

0.27 

0.14 

0.0091 

0.89 

0.47 

0.1014 

4.35 

0.50 

1.9139 

fl 

2.69 

1.04 

0.8321 

2.45 

0.69 

0.6496 

9.73 

0.65 

9.5151 

2.30 

0.59 

0.5612 

4.38 

10.15 

12.225 

9.22 

0.37 

8.5214 

10 

0.15 

1.13 

0.1288 

14.03 

21.95 

67.882 

1.93 

2.44 

0.9664 

11 

3.36 

4.98 

3.6090 

22.54 

9.89 

60.572 

9.50 

9.54 

18.118 

12 

0.00 

2.47 

0.6082 

18.51 

17.95 

66.497 

20.26 

0.44 

41.064 

Table  4.  Prediction  of  soil  parameters  using  the  orthogonal  response  surface  developed  when  a  third  order  correction 

was  employed 
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Table  5.  Sinkage  parameters  obtained  from  filed  tests  conducted  in  a  Yolo  loam  soil. 


Test 

Date 

Soil* 

Condition 

Plate 

Size,  mm 

Sinkage** 

Constants 

kPa 

Sinkage 

Constant,  n 

Overall 

R2 

November, 

Undisturbed 

50.8 

695.6 

0.609 

0.774 

1991 

C=  1 1 .5  kPa 

76.2 

819.6 

0.84 

0.921 

<J>=32.9  deg. 

101.6 

1091.2 

0.828 

0.911 

MC=8.9% 

127.0 

778.3 

0.767 

0.743 

p=1510,  kg/m^ 

e=0.755 

152.4 

1054.4 

0.971 

0.942 

September, 

Undisturbed 

50.8 

959.5 

0.646 

0.776 

1992 

C=32.3  kPa 

$=  212  deg. 

76.2 

921 

0.499 

0.715 

MC  =  5.14% 

p=1510,  kg/m^ 

e=0.755 

101.6 

1839.4 

0.992 

0.934 

Tilled 

50.8 

607.6 

0.918 

0.776 

C=22.7  kPa 

<t>=22.8  deg. 

76.2 

605.6 

0.776 

0.875 

MC=4.55% 

p=1433,  kg/m^ 

e=0.85 

101.6 

796.7 

0.958 

0.875 

*  C  =  cohesion;  4>  =  soil  internal  angle  of  friction;  MC  =  moisture  content,  dry  basis;  p=  bulk  density; 
e  =  void  ratio. 

**  Logarithmic  mean  of  all  eight  replicates. 
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Test  Description 

Best  Results 

Reasonable  Results 

No.  1 

Test  #1 

k=0.0545 

k=0.0446 

Undisturbed  Soil 

*0=0.1310 

"0=0.2089 

November,  1991 

K=0.7124 

K=0.8928 

C=  12.348  kPa 

C=1 1.371  kPa 

4>=22.5  deg. 

<>=22.5  deg. 

SSE=0.00380 

SSE=0.00490 

SSR=0.00249 

SSR=0.00601 

Test  #2 

k=0.0513 

k=0.0225 

Undisturbed  Soil 

*0=0.159 

b=0.370 

September,  1992 

K=0.7365 

K=0.8950 

0=14.321  kPa 

C=15.068  kPa 

<>=28.333  deg. 

<>=26.050  deg. 

SSE=0.00665 

SSE=0.01332  * 

SSR=0.00089 

SSR=0.00345 

Test  #3 

k=0.1037 

k=0.0505 

Tilled  Soil 

*0=0.1859 

*0=0.370 

September  1992 

K=0.6247 

K=0.8993 

C=8.997  kPa 

C=16.516  kPa 

<>=31.820  deg. 

<>=25.247  deg. 

SSE=0.00378 

SSE=0. 00444 

SSR=0.00016 

SSR=0.00337 

Reasonable  Results 
No.  2 


k=0.0638 
d=0.0983 
K=0.7955 
C=15.344  kPa 
<>=22.5  deg. 
SSE=0.02146  * 
SSR=0.01295 


*  Relatively  high  SSE  indicating  an  unreasonable  solution. 
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Figure  4.  The  plot  of  error  without  any  correction  in  Eq  -  Ej  space 


Si  T  w-  go 


Figure  5.  The  plot  of  orthogonal  response  surface  in  Eq  -  Ej  space  when  a  second  order 

correction  was  included 


Figure  7.  Plot  of  orthogonal  response  surface  in  Eq-  Et  space  when  a  third  order 

correction  was  included 
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Figure  8.  Plot  of  error  in  Eq  -  Et  space  when  a  third  order  correction  was  included 


Nondimensional  Plate  Sinkage 
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Nondimensional  Parameter  Values 


Figure  9.  Plate  sinkage  vs.  nondimensional  values  of  parameter  k  for  a  4  in  plate  subjected 

to  80  PSI  applied  load 


Nondimensional  Plate  Sinkage 
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R2  =0.9927 


Nondimensional  Parameter  Values 

Figure  11.  Plate  sinkage  vs.  nondimensional  values  of  parameter  K  for  a  4  in  plate 

subjected  to  80  PSI  applied  load 
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R2  =1.0000 


Nondimensional  Parameter  Values 

Figure  12.  Plate  sinkage  vs.  nondimensional  values  of  parameter  C  for  a  4  in  plate 

subjected  to  80  PSI  applied  load 
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R2  =1.0000 


Nondimensional  Parameter  Values 

Figure  13.  Plate  sinkage  vs.  nondimensional  values  of  parameter  4>  for  a  4  in  plate 

subjected  to  80  PSI  applied  load 
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R2  =1.0000 


Nondimensional  Parameter  Values 

Figure  14.  Plate  sinkage  vs.  nondimensional  values  of  parameter  e  for  a  4  in  plate  subjected 

to  80  PSI  applied  load 


Response  Surface 
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Real  Surface 


Figure  15.  Response  surface  vs.  real  surface  for  a  4  in  plate  subjected  to  15  PSI  pressure  - 

without  correction 


y  =  1.009x  -  0.013  ^  =  0.997 


Real  Surface 


Figure  16.  Response  surface  vs.  real  surface  for  a  4  in  plate  subjected  to  SO  PSI  pressure 

without  correction 


47 


Figure  17.  Response  surface  vs.  real  surface  for  a  4  in  plate  subjected  to  15  PSI  pressure  - 

using  second  order  correction 
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y  =  1.018x  -  0.023  r2  =  0.997 


Figure  18.  Response  surface  vs.  real  surface  for  a  4  in  plate  subjected  to  80  PSI  pressure  - 

using  second  order  correction 
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Figure  19.  Response  surface  vs.  real  surface  for  a  4  in  plate  subjected  to  15  PSI  pressure  - 

using  third  order  correction 
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Figure  20.  Response  surface  vs.  real  surface  for  a  4  in  plate  subjected  to  80  PSI  pressure  - 

using  third  order  correction 
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y  =  -  0.37058  +  2.2590x  RA2  =  0.989 


Figure  23.  Experimental  and  simulated  sinkage  for  a  2  in  plate  in  an  undisturbed  soil 
when  the  soil  parameters  obtained  using  the  inverse  solution  technique  using  a  four  inch 
plate  test  was  used  to  simulate  2  in  plate  behavior.  Experimental  results  correspond  to 
September  1992  tests. 
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SOIL  STRESSES  UNDER  TRACTOR  TIRES 


A.C.  Bailey,  R.L.  Raper,  C.E.  Johnson,  T.R.  Way,  and  E.C.  Burt.1 

Introduction 

Soil  stress  state  transducers  (SST's)  have  been  used  to  determine  the  state  of  stress 
beneath  tractor  tires  operating  under  a  range  of  dynamic  loads  and  inflation  pressures.  The 
stress  state  transducer  has  been  described  by  Nichols  et  ai.(1987).  These  SST's  provide  the 
data  necessary  to  calculate  the  complete  stress  state  at  the  transducer.  The  stress  state  is  a  3 
x  3  symmetric  tensor  with  6  independent  stresses,  and  may  be  represented  by  the  6 
independent  stresses,  the  three  principal  stresses  (o„  o2,  o3)  and  their  directions,  or  stresses 
on  particular  planes,  such  as  the  octahedral  shearing  (t**)  and  normal  (a,**)  stresses  and  their 
directions  (Bailey  and  Burt  1988). 

Procedure 

A  recent  experiment  (Bailey  et  al.  1993)  studied  the  effect  on  soil  stresses  of  4 
combinations  of  tire  loads  and  inflation  pressures,  Table  1 .  Treatments  L  and  H  are 
combinations  of  load  and  inflation  pressures  taken  from  the  manufacturer's  recommendations 
(Goodyear  1992).  Treatment  U  (underload,  load  less  than  recommendation  for  the  inflation 
pressure)  used  the  load  from  the  L  level  (13.1  kN)  and  the  inflation  pressure  from  the  H 
treatment  (124  kPa).  Treatment  O  (overload)  treatment  used  the  load  from  the  H  level  (25.3 
kN)  and  the  inflation  pressure  of  the  L  treatment  (41.4  kPa).  The  tire  was  operated  at  a 
constant  forward  velocity  of  0.15  m/s  and  a  constant  slip  of  10%. 

Two  soils,  Norfolk  sandy  loam  (NSL)  and  Decatur  clay  loam  (DCL),  and  two  profiles 
in  each  soil  were  used.  One  profile  was  relatively  loose  and  uniform.  A  hardpan  was 
present  in  the  second  profile.  A  SST  was  placed  at  the  depth  of  the  hardpan  (hardpan  depth) 
and  a  second  SST  placed  midway  between  the  surface  and  the  hardpan  (shallow  depth).  Both 
SST's  were  aligned  directly  under  the  centerline  of  the  path  of  the  tire.  Soil  bulk  density 
samples  were  taken  in  each  tireprint  at  the  depth  of  each  SST  after  completion  of  the  tests. 

Results  and  Discussion 

Typical  data  from  a  SST  are  shown  in  figure  i.  The  data  shown  are  from  the  H 
treatment  at  the  hardpan  depth  in  the  NSL  with  a  hardpan.  The  value  of  0  on  the  distance 
axis  represents  the  horizontal  location  of  the  tire  axle,  and  data  at  positive  distances  represent 


‘The  authors  are:  Alvin  C.  Bailey  and  Randy  Raper,  Agricultural  Engineers,  National 
Soil  Dynamics  Laboratory,  USDA-ARS,  Auburn,  AL.;  Clarence  E.  Johnson,  Professor, 
Agricultural  Engineering  Dept.,  Alabama  Agricultural  Experiment  Station,  Auburn 
University,  AL.;  and  Thomas  R.  Way,  Agricultural  Engineer,  and  Eddie  C.  Burt,  Research 
leader,  National  Soil  Dynamics  Laboratory,  USDA-ARS,  Auburn,  AL. 
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Figure  1 .  Typical  measured  pressures  from 
the  stress  state  transducer  under  the 
centerline  of  the  tire. 


Distance,  m 


Figure  2.  Calculated  stresses  from  data 
shown  in  figure  1. 


pressures  in  front  of  the  tire.  Figure  2  presents  the  calculated  principal  stresses  (a,,  a2,  and 
a3)  and  the  octahedral  stresses  (a**  and  x^  for  the  data  from  figure  1  using  the  same 
distance  axis.  The  major  principal  stress,  a„  is  the  greatest  stress  but  displays  the  same 
trends  as  the  octahedral  shearing  stress.  For  this  discussion  the  octahedral  stresses  and 
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will  be  used  to  represent  the  stress  state  in  the  soil  under  the  tractor  tire.  The  peak  values 
of  Con  and  the  corresponding  values  of  were  selected  for  further  analyses.  Table  1 
presents  of  the  means  of  the  4  replications  of  and  x*,  from  all  treatments,  soils,  and 
profiles. 

The  peak  values  of  and  the  corresponding  values  of  were  analyzed  with  SAS 
(1990)  using  a  factorial  design  on  each  soil  separately.  Replications  were  nested  within  soil 
profile.  The  multivariate  option  MANOVA  was  also  used  to  analyze  both  octahedral  stresses 
together. 

Higher  levels  of  either  dynamic  load  or  inflation  pressure  generated  higher 
octahedral  stresses  beneath  the  tractor  tire  in  both  soils  when  averaged  across  all  depths  and 
soil  profiles  (figures  3  and  4).  Both  dynamic  load  and  inflation  pressure  were  significant 


factors  (5%  level)  affecting  the  two  octahedral  stresses,  except  for  inflation  pressure  on  in 
the  DCL.  This  exception  is  probably  because  the  DCL  had  the  greatest  variability  in  x**  (a 
deviatoric  stress  contrasted  to  a**,  an  average  stress).  Both  dynamic  load  and  inflation 
pressure  were  significant  factors  (5  %  level)  affecting  the  bulk  density  and  increase  in  bulk 
density  at  each  depth  level  in  each  soil. 

The  mean  net  traction  and  tractive  efficiency  data  for  all  treatments  in  each  soil  are 
presented  in  Table  2.  An  ANOVA  of  the  net  traction  and  tractive  efficiency  showed  that 
inflation  pressure  and  dynamic  load  both  had  significant  effects  on  both  performance 
variables,  and  that  higher  inflation  pressures  had  lower  net  traction  and  tractive  efficiencies  at 
the  same  dynamic  load.  These  results  support  the  conclusions  from  the  soil  stresses  and  bulk 
density  data.  Higher  inflation  pressures  at  the  same  dynamic  load  have  lower  tractive 
efficiencies  and  generate  higher  soil  stresses  and  soil  compaction. 


Table  2.  Mean  tractive  performance  data. 
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Soil  Compaction  Research  Needs 

P.  T.  Corcoran1 


Introduction 

There  are  two  needs  for  compaction  research  that  are  potential  topics  for  the  workshop  on 
"Modeling  the  Mechanics  of  Off-Road  Mobility".  Those  topics  are  predictive  analytical  models 
for  machine  compaction,  and  the  forecasting  of  machine  compaction  performance  on  specific  job 
site  conditions. 

These  two  needs  are  similar  in  that  both  are  needed  to  reduce  the  requirements  for  machine 
testing.  Analytical  models  for  machine  compaction  can  reduce  the  need  for  expensive  prototype 
testing,  and  improve  the  understanding  of  relationships  between  machine  and  soil  in  the 
production  of  compaction.  Performance  forecasting  can  eliminate  the  need  for  proof  testing,  and 
improve  the  quality  of  soil  structures. 


Analytical  Compaction  Models 

There  has  long  been  a  need  for  accurate  and  complete  analytical  compaction  models.  The  ideal 
is  to  have  such  models  based  on  proven  theory,  and  not  be  dependent  on  relationships  defined 
only  by  empirical  testing.  Although  empirical  testing  is  often  the  most  attractive  as  a  low  risk 
practical  means  to  define  machine/soil  relationships,  empirical  testing  is  always  inherently 
limited  in  application  to  those  parameters  and  ranges  of  parametric  values  included  in  whatever 
experiments  provide  the  data  base  for  empirical  relationships. 

Conversely,  theoretically  based  analytical  models  are  essentially  unlimited  in  range  of 
application.  The  more  complete  the  theory  and  the  more  detailed  the  model,  the  broader  the 
range  of  application.  Historically  analytical  soil  compaction  models  have  been  limited  by  a  lack 
of  computational  capacity  and  detail  to  handle  the  non-linear  and  fine  grained  characteristics  of 
machine  soil  compaction  problems.  These  technological  hurdles  are  coming  down  with  the 
computational  speeds  of  modem  computer  work  stations,  and  the  sophistication  in  newer 
modeling  software. 

Recent  work  with  a  finite  element  soil  compaction  model  is  confirming  the  technological  hurdles 
to  a  comparatively  complete  and  detailed  model  are  gone  or  greatly  reduced.  The  challenge  now 
becomes  one  of  using  this  technological  opportunity.  The  author  is  interested  in  sharing 
experiences  with  others  addressing  the  need  for  predictive  models,  both  analytical  and  empirical. 
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Compaction  Performance  Forecasting 

Compaction  is  very  often  a  critical  need  in  the  construction  industry.  Achieving  compaction 
specifications  can  be  the  critical  path  guiding  the  progress  of  a  construction  site,  and  therefore 
achieving  compaction  specification  rrr  the  determinant  factor  in  scheduling  and  costing  a 
job.  Contractors  need  accurate  inforn  on  the  ability  to  meet  compaction  specifications,  and 
traditionally  this  has  required  proof  tt  ,.g.  Proof  testing  can  be  far  from  an  accurate  means  to 
insure  compaction  specifications  will  actually  be  met.  The  common  variability  in  job  site  soils, 
virgin  soil  conditions,  and  the  influence  of  weather  can  destroy  the  accuracy  of  proof  testing. 

An  ideal  would  be  the  ability  to  forecast  the  performance  of  compaction  machines  based  on 
specific  job  site  conditions,  and  adapt  machines  and  operating  procedures  to  the  job  site  as 
needed  for  both  maximum  efficiency  and  performance.  However,  most  contractors  lack 
information  and  interpretation  to  achieve  such  an  ideal. 

There  appear  to  be  two  major  needs  to  allow  for  the  iri’r  "'ements  of  performance  forecasting 
and  the  reduction  of  need  for  proof  testing;  the  quanL  .t  .  .  onships  between  machine 
parameters,  job  site  conditions,  and  compaction  specifications;  and  the  relationship  between 
laboratory  compaction  measurement  and  actual  field  capability. 

Quantification  of  relationships  between  machine,  soil,  and  specifications  is  perhaps  little  more 
difficult  then  recording  information  already  available  from  most  job  siies  to  establish  a  data  base. 
There  may  be  some  incremental  additional  information  required,  however,  it  would  be  contrary 
to  the  purpose  of  performance  forecasting  to  require  significant  additional  information  than  "hat 
normally  available.  A  requirement  for  significant  additional  information  would  only  trade-off 
one  inefficiency  for  another  and  thus  not  bring  about  a  net  gain  in  overall  job  efficiency.  A 
reasonably  inclusive  database  should  then  give  the  opportunity  to  establish  correlations  yielding 
a  forecasting  capability. 

Improving  the  relationship  between  laboratory  compaction  measurement  and  actual  field 
capability  may  be  more  difficult.  The  problem  is  compaction  machines  in  the  field  do  not  have 
the  same  energy  efficiency  nor  use  the  same  mechanism  to  produce  soil  compaction  as  standard 
laboratory  procedures.  Therefore,  densities  achievable  in  the  laboratory  are  not  necessarily 
achievable  in  the  field  at  the  same  energy  level  and  maybe  not  at  any  energy  level  depending  on 
compaction  equipment  available.  Additionally,  as  shown  by  energy/density  relationships  from 
laboratory  testing,  optimum  moisture  levels  vary  with  compaction  energy  and  may  vary  with 
method  of  compaction.  The  need  is  to  provide  information  to  guide  field  use  of  compaction 
machines  both  for  the  selection  of  the  optimum  machine,  and  the  optimum  use  of  the  machine. 
Obtaining  such  information  could  require  a  significant  amount  of  controlled  testing  and  therefore 
be  a  relatively  costly  endeavor. 


One  possible  alternative  to  establishing  correlation  between  laboratory  and  field  through 
empirical  testing  would  be  the  development  and  verification  of  fundamental  theory  of  soil 
compaction.  Such  a  theory  may  already  be  in  existence  based  on  the  compaction  energy /density 
relationships  defined  by  current  laboratory  testing.  Therefore,  expansion  of  a  theory  of  this  type 
to  correlation  between  laboratory  energy  and  actual  field  energy  could  become  the  foundation  for 
performance  forecasting  and  not  require  large  amounts  of  machine  tests. 

There  may  be  a  significant  amount  of  compaction  forecasting  information  already  available  but 
not  widely  disseminated,  or  there  may  be  a  significant  amount  of  information  available  from  a 
wide  distribution  of  sources  but  not  consolidated,  or  there  may  simply  be  a  lack  of  information 
pertinent  to  the  need.  The  author  is  interested  in  sharing  thoughts  and  ideas  on  the  need,  value, 
and  feasibility  for  accurate  performance  forecasting  of  soil  compaction  and  the  availability  of 
pertinent  information. 


Localized  Energy  Dissipation  in  Strained  Granular  Materials 

Peter  K.  Haff» 


A  granular  material  is  a  mechanical  system  composed  of  distinct  macroscopic 
interacting  components.  A  soil,  a  sand  dune,  a  heap  of  mine  tailings,  and  a  fractured  rock 
mass,  are  all  granular  materials,  with  "grain  sizes"  ranging  from  microns  to  meters  in 
diameter.  There  are  two  approaches  that  can  be  used  to  model  granular  systems 
quantitatively.  A  continuum  method  based  upon  a  partial  d&'erential  equation  or  a  discrete 
method  that  retains  specific  reference  to  the  particulate  nature  of  the  medium. 

Continuum  models  require  the  existence  of  an  "averaging  volume",  large  compared 
to  grain  sizes,  within  which  the  characteristics  of  individual  particles  can  be  replaced  by 
variables  like  velocity,  density,  components  of  stress,  etc.  Discrete  models  on  the  other 
hand  retain  reference  to  individual  particle  identities.  In  discrete  models  one  can  investigate 
"microscopic"  behavior  that  is  lost  or  obscured  in  the  averaging  transition  to  a  continuum 
model.  One  such  type  of  behavior  in  granular  systems  is  the  nature  of  frictional  losses 
incurred  when  the  medium  is  subjected  to  inelastic  strain. 

A  drawback  in  using  a  discrete  modeling  technique  is  that  the  number  of  discrete 
particles  that  can  be  handled  is  limited  (thousands  to  tens  of  thousands  of  particles),  and 
hence  we  can  study  only  a  small  volume  of  material  at  one  time,  while  with  a  continuum 
model  we  can  usually  model  large  volumes.  Consequently,  discrete  techniques  are  often 
best  used  to  generate  insight  into  microscopic  mechanisms  that  can  help  us  in  interpretation 
of  larger  scale  modeling,  or  to  help  construct  constitutive  relations  for  use  in  continuum 
models. 

In  a  granular  medium  subjected  to  a  given  load  (a  generic  prototype  of  a  soil  or 
similar  granular  material  subjected  to  vehicle  loading),  the  stress  is  distributed  throughout 
the  medium  via  grain  contacts  and  fluid  pressure  forces.  For  simplicity,  we  consider  here 
only  dry,  noncohesive  materials.  As  strain  develops  in  the  material,  compression  of 
individual  grains  occurs  at  the  grain-grain  contacts.  If  the  contacts  do  not  slip,  and  the 
strain  rate  is  small,  then  the  strain  is  generally  reversible  upon  unloading  and  energy  loss  is 
zero.  At  larger  strains,  grain  contacts  can  slip.  Since  the  contacting  surfaces  of  earth 
materials  are  frictional.  Coulomb-type  losses  are  incurred,  and  contact  positions  do  not 
return  to  their  original  configurations  upon  unloading. 

Particle  dynamics  studies  of  energy  dissipation  in  strained  granular  material  point  to 
the  importance  of  fluctuations  in  the  stress  distribution  at  grain  contacts.  These  fluctuations 
can  influence  the  macroscopic  lossiness  of  the  material.  When  a  compressive  load  is 
applied  to  a  granular  assembly,  the  allocation  of  stress  among  contacts  is  determined  in 
large  measure  by  the  geometric  placement  of  grains.  Consider  for  simplicity  a  system  of 
circular  or  spherical  particles.  Each  particle  has  a  set  of  neighboring  particles  with  whom  it 
is  in  contact  By  drawing  an  imaginary  line  between  the  centers  of  every  pair  of  contacting 
grains,  we  define  a  stress  network.  A  normal  force  associated  with  grain  stiffness  is 
associated  with  each  element  in  this  network,  as  well  as  a  tangential  force  due  to  friction 
between  the  contacting  grains.  Die  normal  force  exerted  across  each  element  of  the  network 
is  a  function  of  the  deformation  at  that  contact  If  the  elastic  interaction  is  modeled  by  a  stiff 
spring,  then  the  force  is  a  function  of  the  instantaneous  compression  of  the  spring.  In 
general  the  stresses  will  be  different  at  each  contact  i.e.,  spring  compressions  will  differ, 
so  that  there  is  a  distribution  of  local  stress  determined  by  the  details  of  grain  packing. 
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Small  amounts  of  strain  can  often  be  accommodated  by  reversible  compression  of 
the  elastic  springs  at  the  contacts,  but  large  strain  must  be  accompanied  by  slippage  and 
rotation  of  grains.  Wherever  slippage  takes  place,  frictional  forces  will  lead  to  energy  loss. 
Because  there  is  a  distribution  of  normal  contact  forces,  some  contacts  will  slip  more  easily 
than  others.  If  a  contact  is  tightly  compressed,  so  that  the  normal  force  is  large,  then  a 
correspondingly  large  local  tangential  stress  is  needed  at  that  contact  to  cause  it  to  fail. 
Conversely,  at  a  contact  where  the  normal  spring  is  not  much  compressed,  slippage  is 
relatively  easy.  The  frictional  energy  lost  in  a  slip  event  is  the  product  of  the  tangential  force 
acting  during  the  slip  and  the  amount  of  displacement  incurred  at  the  contact  Assuming  for 
the  sake  of  argument  equal  displacements,  die  greatest  amount  of  energy  lost  per  slipping 
contact  is  at  the  strong  contacts.  Conversely,  the  smallest  amount  of  energy  lost  per  slip 
event  is  at  the  weakly  compressed  contacts. 

Particle-dynamics-model  simulations  of  sheared  granular  materials  suggest  that  for 
the  system  as  a  whole,  most  energy  is  lost  at  contacts  of  intermediate  strength.  The  strong 
contacts  do  not  slip  sufficiently  often  to  dominate  the  overall  energy  loss,  while  the  weak 
contacts  slip  often  but  do  not  generate  enough  loss  per  contact  to  be  the  dominant  loss 
mechanism.  A  consequence  of  this  observation  is  that  the  macroscopic  rate  of  energy  loss 
in  deformation  is  not  a  simple  function  of  the  grain-on-grain  frictional  properties  of 
individual  grains.  The  energy  loss  will  usually  be  less  than  that  expected  from  a  simple 
averaging  approach,  the  magnitude  of  the  effect  depending  upon  the  detailed  distribution  of 
forces  over  the  stress  network.  A  corollary  is  that  since  strong  contacts  cannot  slip  easily, 
clusters  of  grains  with  strong  contacts  between  them  tend  to  rotate  as  a  whole,  with  slip  and 
energy  loss  occurring  on  the  periphery  of  the  cluster.  As  the  material  strains,  the  stress 
network  continually  adjusts  to  maintain  force  balance.  Rotating  friction-locked  clusters 
eventually  unlock  as  their  internal  stress-network  becomes  less  well  aligned  with  the  overall 
stress  field  in  the  medium.  At  this  point  the  local  stresses  readjust,  slippage  begins  to  occur 
within  the  previously  rigid  cluster,  and  new  clusters  spontaneously  appear  nearby. 

The  frictional  losses  and  hence  the  energy  absorption  of  the  macroscopic  granular 
medium  depend  upon  the  details  of  such  microscopic  mechanisms.  By  elucidating  these 
mechanisms,  discrete  computational  models  can  help  to  further  our  understanding  of  the 
dynamic  response  of  granular  systems  such  as  soils  to  external  loads. 
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A  Case  for  Improved  Soil  Models  in  Tracked  Machine  Simulations 

F.B.  Huck' 


Abstract 

A  planar,  multibody  dynamics  model  of  a  track-type-tractor  was  developed  to  provide  an 
analytical  tool  to  evaluate  alternative  new  tractor  designs  and  to  resolve  problems  on  current 
products.  The  model  has  been  applied  to  studies  of  track  chain  vibrations,  track/sprocket 
jumping,  sprocket/bushing  wear,  and  the  influence  of  track  and  undercarriage  kinematics  on 
fundamental  tractor  rigid  body  vibration  mode  excitation.  Applications  of  the  model  in  studies 
of  higher  frequency  vehicle  vibrations  associated  with  ride  quality  demonstrated  the  important 
influence  that  track/soil  interface  models  have  on  machine  vibration  excitation  and  the  need  to 
improve  our  capabilities  in  this  area. 

Introduction 

Earthmoving  equipment  manufacturers  must  identify  ways  to  reduce  product  design  and 
development  time  and  cost  to  remain  competitive  in  world  markets.  Manufacturers  increasingly 
turn  to  engineering  computer  analysis  as  a  time  and  cost  savings  methodology  to  supplement  and 
ultimately  reduce  dependence  on  their  more  traditional  and  costly  build  and  test  approach  to 
product  design  and  development. 

As  computers  continue  to  increase  in  speed  and  expand  in  memory  capacity,  engineering 
software  functionality  grows  to  quickly  fill  any  vacant  memory  cell  or  unused  CPU  cycle. 
Performance  analysts,  as  a  result,  are  able  to  develop  increasingly  realistic  and  detailed 
interdisciplinary  models  to  simulate  the  overall  dynamic  response  of  the  complete  earthmoving 
machine. 

The  response  of  the  machine  ultimately  depends  upon  the  external  forces  which  act  upon 
it.  In  the  case  of  an  earthmoving  machine,  these  include  gravitational  forces,  combustion  of 
engine  fuel,  operator  interactions,  and  soil  reactions.  None  of  these  external  influences,  with  the 
exception  of  gravity,  is  as  well  understood  from  an  analytical  standpoint  as  it  needs  to  be  to 
match  the  degree  of  sophistication  now  attainable  in  dynamics  models  of  the  machine  itself.  Of 
the  three  least  understood  forces,  the  soil  force  reaction,  though  possibly  not  as  analytically 
intractable  as  the  human  operator,  is  the  most  critically  lacking  component  in  most  earthmoving 
machine  performance  models. 
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Discussion 


The  Model 

One  example  of  a  model  that  pushed  the  limits  of  computing  resources  and  analysis 
software  to  their  limits  at  the  time  of  its  development  was  a  2-  dimensional,  pitch  plane,  multibody 
dynamics  model  of  a  high-drive  track-type-tractor.  Initiated  in  the  mid  1980's,  the  model  was 
developed  within  the  framework  of  the  commercial,  rigid  multibody  dynamics  code,  DRAM 
(Chase  and  Angell  1977).  The  model  (Figures  1  &  2)  (Huck,  1987)  was  unique  from  previous 
tracked  machine  models  in  that  it  treated  each  link  in  the  track  chain  as  a  distinct  rigid  body.  This 
level  of  model  fidelity  was  required  to  address  the  range  of  questions  being  asked  at  the  time. 

Typical  Model  Applications 

The  model  was  applied  primarily  to  questions  related  to  track  chain  and  undercarriage 
dynamics  or  to  the  influence  of  track  and  undercarriage  dynamics  on  fundamental  machine 
vibration  modes  excitation.  In  the  area  of  track  chain  dynamics,  the  model  was  used  to 
understand  the  relationship  among  track  pitch,  catenary  length,  track  tension,  and  track  speed  on 
the  excitation  of  transverse  vibration  modes  in  the  track  catenaries. 

A  second  application  demonstrated  how  premature  wear  develops  when  out-of-tolerance 
sprocket  teeth  segments  lead  to  adverse  track  bushing  and  sprocket  tooth  engagement.  A  third 
application  simulated  the  track/sprocket  tooth  jumping  that  can  occur  during  rapid 
forward/reverse  directional  shifts. 

The  primary  application  for  the  model  is  to  optimize  undercarriage  component  placement 
in  order  to  minimize  the  influence  of  track  and  undercarriage  kinematics  on  the  excitiation  of  the 
chassis'  fundamental  rigid  body  pitch  vibration  mode  in  the  4  to  6  Hz  range.  If  overly  excited, 
this  mode  can  affect  dozer  controllability  during  finish  grading  applications. 

Importance  of  Soil  Models 

The  DRAM  tracked  machine  model  incorporates  simple  representations  of  soil  behavior 
to  simulate  the  ground  reaction  forces  that  support  and  propel  the  machine. 

Early  versions  of  these  models  applied  Bekker  pressure-sinkage  like  relationships  to 
simulate  the  normal  support  force  acting  on  each  individual  track  shoe.  (Bekker,  1969)  Likewise, 
a  model  similar  to  that  proposed  by  Kacigen  &  Guskov  was  applied  to  each  track  shoe  grouser  to 
simulate  the  soil  shearing  forces  associated  with  traction.(Kacigen  &  Guskov,  1968)  These 
models  are  simply  nonlinear  elastic  representations  of  the  soil.  There  is  no  coupling  between  the 
normal  support  and  traction  force  models  other  than  the  dependence  of  the  tractive  force 
generated  by  a  given  track  shoe  on  the  normal  pressure  imposed  on  the  soil  by  the  adjacent  shoe. 

While  such  models  may  be  adequate  for  steady  state  predictions  of  tracked  machine 
pull-slip  characteristics,  they  are  quite  inadequate  for  detailed  dynamic  models  of  the  kind 
discussed  here.  Certainly  a  major  deficiency  is  the  inability  for  a  nonlinear  elastic,  single 


modulus,  spring-only  model  to  provide  the  energy  dissipation  in  the  soil  that  accompanies  the 
passage  of  a  tracked  machine. 

As  mentioned  above,  the  tracked  machine  model  is  commonly  used  to  predict  machine 
vibrations.  Numerous  studies  over  the  years  show,  not  unexpectedly,  that  its  predictive  accuracy 
is  strongly  dependent  upon  how  the  soil  is  characterized. 

The  plots  in  Figure  3  show  the  horizontal  component  of  machine  velocity  predicted  for  the 
CG  of  a  D8L  tractor  operating  at  a  sprocket  speed  of  1 .89  MPH  (0.845  m/s)  under  high  drawbar 
load.  The  two  simulations  are  identical  except  for  the  presence  or  absence  of  an  ad  hoc  damping 
term  that  was  added  to  the  traction  force  model.  The  top  plot,  from  the  simulation  with  no 
traction  damping,  shows  the  superimposition  of  a  1  Hz  and  a  4  Hz  vibration  in  forward  tractor 
motion.  The  1  Hz  component  is  the  result  of  a  1  Hz  variation  in  track  tension.  It  is  consistent 
with  the  natural  frequency  of  the  machine  mass  "bouncing"  horizontally  against  the  nonlinear 
traction  "springs".  The  4  Hz  component  (actually  3. 8 8  Hz)  represents  the  track  first  order 
chordal  excitation.  The  introduction  of  a  50%  critical  damper  into  the  traction  model  resulted  in 
the  response  shown  on  Figure  3b.  The  vibration  content  of  the  machine  is  now  almost  purely  due 
to  chordal  excitation  as  would  be  expected. 

Similar  results  are  obtained  in  predictions  of  vertical  machine  motion.  As  a  result  of  these 
early  model  development  observations,  ad  hoc  viscous  damping  terms  were  added  to  both  the 
normal  and  tractive  force  components  of  the  track/soil  interface  model  to  account  for  energy 
absorption  by  the  soil. 

Soil  Model  Influence  on  Machine  Vibration  Predictions 

During  the  course  of  normal  use,  track  link  rails  will  wear  and  develop  a  scalloped  profile. 
The  ride  quality  of  the  machine  can  be  adversely  effected  if  the  amplitude  of  the  scallop  pattern 
becomes  excessive.  The  pattern  is  usually  adequately  described  as  a  Fourier  series  sine  wave 
expansion  to  third  order  of  track  pitch.  This  pattern  is  well  known  to  undercarriage  designers; 
and,  they  normally  account  for  it  by  properly  positioning  the  track  rollers. 

In  rare  instances  on  certain  soils,  scallop  patterns  with  fourth  order  content  appear.  A 
recent  case  provided  the  opportunity  not  only  to  apply  the  DRAM  tracked  machine  model  to  an 
analysis  of  the  situation  but  also  to  calibrate/validate  the  model. 

Full  scale  machine  vibration  tests  on  the  subject  factor  were  run  to  accumulate  data  on 
the  fourth  order  vibration  phenomenon.  Triaxial  accelerometers  were  positioned  at  the  front  and 
rear  of  the  roller  frames  (RF),  and  at  the  front  and  rear  of  the  main  frame  (MF)  or  chassis.  Fast 
Fourier  Transform  (FFT)  plots  of  the  measured  vertical  acceleration  data  clearly  show  the  fourth 
order  content  in  the  signal,  particularly  on  the  roller  frame  signals.  (See  Figure  4) 

A  DRAM  tracked  machine  model  was  assembled  for  the  subject  tractor  an-  simulations 
run  to  predict  accelerations  at  points  where  they  were  measured  on  the  machine.  Frequency 
response  plots  of  the  predicted  accelerations  are  shown  in  Figure  5.  Note  that  the  only  model 
parameters  changed  for  the  six  simulation  results  shown  were  the  effective  spring  and  damping 
coefficients  in  the  normal  soil  support  force  model.  Clearly,  soil  parameter  changes  have  a  strong 
influence,  not  only  on  the  relative  magnitude  of  a  particular  order  but  also  on  the  order  that 
predominates. 


The  simulation  results  with  a  Bekker  coefficient  of  2.0E6  and  a  damping  coefficient  of 
6.0E4  most  closely  resemble  the  measured  data.  These  parameters  were  chosen  to  "calibrate"  the 
soil  model  for  subsequent  work  on  this  project.  Fourth  order  scalloped  wear  patterns  are  now 
routinely  considered  during  undercarriage  design  studies  that  employ  this  model  to  determine 
track  roller  placement 

Improved  Transient  Soil  Models 

The  work  discussed  above  is  just  one  example  of  the  need  to  improve  analytical 
descriptions  of  the  soil's  transient  response  characteristics.  One  would  like  more  realistic  models 
that  can  be  simply  and  independently  calibrated.  It  is  far  too  costly  and  time  consuming  to 
conduct  full  scale  tests,  like  the  one  above,  to  calibrate  the  model.  Furthermore,  it  defeats  a  main 
purpose  of  analysis,  which  is  to  provide  a  reliable,  predictive  capability  without  the  need  to  build 
and  test  a  prototype. 

To  this  end.  Caterpillar  has  worked  to  develop  improved  track  soil  interface  models.  A 
new,  semi-empirical  visco-elasto-plastic  model  that  more  accurately  characterizes  the  plastic 
deformation  that  occurs  during  dynamic,  repetitive  soil  loading  has  been  developed  by  Hombrook 
(Figure  6)  (Hombrook  1992).  This  "dual  stiffness"  spring  and  viscous  damper  model  is  simple  to 
implement,  relatively  easy  to  calibr  ,te,  and  does  a  very  good  job  of  matching  measured 
force/sinkage  characteristics  of  repetitively  loaded  flat  plates  in  plastic  soil  with  frictional  and 
cohesive  characteristics.  (Figure  7) 

The  new,  visco-elasto-plastic  model  has  been  implemented  in  the  DRAM  tracked  machine 
model;  but,  it's  impact  on  the  total  model's  predictive  accuracy  in  machine  vibration  simulation 
studies  is  yet  to  be  verified. 

Similar  work  is  needed  to  improve  the  tractive  force  portion  of  the  track/soil  interaction 
model.  To  date,  little  effort  has  been  devoted  to  that  task. 

Summary 

Earthmoving  machine  manufacturers  increasingly  depend  upon  sophisticated  engineering 
analysis  techniques  to  help  reduce  the  time  and  costs  involved  in  new  product  introduction. 
Advances  in  computer  technology  now  permit  the  development  of  detailed  machine  performance 
models  with  the  capability  to  address  a  broad  range  of  design  and  development  questions.  Model 
accuracy  must  approach  test  measurement  accuracy  if  costly  and  time  consuming  full  scale 
testing  is  to  be  minimized.  However,  earthmoving  machine  performance  model  accuracy  is 
strongly  dependent  on  methods  employed  to  characterize  machine/soil  interactions.  While  some 
progress  has  been  made,  there  remains  an  urgent  need  for  computationally  efficient,  easy  to 
calibrate  models  which  accurately  describe  the  transient  response  characteristics  of  soil. 
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Model  Predictions  of  Scalloped  Track  Induced  Frame  Accelerations 
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Prediction  of  Soil  Compaction  Behavior 

Clarence  E.  Johnson,  Alvin  C.  Bailey,  and  Randy  L.  Raper' 


Introduction 

Our  research  interests  include  the  development  of  mathematical  models  of  soil  compaction 
behavior  (constitutive  relations)  of  agricultural  topsoil  necessary  for  prediction  of  stress  propagation  in 
the  soil  and  resulting  state  of  compactness.  The  source  of  the  force  systems  which  influence  compaction 
may  come  from  field  machinery  moving  over  the  soil  surface  or  through  the  soil  beneath  die  surface. 

The  initial  state  of  die  soil  may  be  very  loose  (following  a  tillage  operation)  and  it  is  often  unsaturated 
either  before  or  after  application  of  the  force  system.  The  state  of  soil  compactness  influences  vegetative 
response,  soil  erosiveness  and  degradation,  and  other  utility  of  our  natural  resources. 

Prevost  (1987)  recognized  that  modern  tools,  such  as  computer  technology  coupled  with  finite 
element  techniques,  provide  die  potential  for  solving  problems  associated  with  soil  behavior  of  greater 
complexity  than  did  past  historical  technology.  But,  he  emphasized  that, 

"Further  progress  in  expanding  analytical  capabilities  in  geomechanics  now  depends 
upon  consistent  mathematical  formulations  of  generally  valid  and  realistic  material 
constitutive  relations." 

The  goal  of  our  research  is  to  develop  a  useful  model  of  soil  stress-strain-strength  behavior  to  predict 
satisfactorily  soil  (and  machine)  "performance"  in  circumstances  important  to  production  agriculture, 
forestry  and  off-road  mobility. 

Current  Status 

Schafer  et  al.  (1989)  summarized  the  status  of  our  soil  compaction  modeling  effort.  Our  current 
deviatoric  stress  model,  described  by  Bailey  and  Johnson  (1989)  is  a  modification  of  a  previous  model 
(Bailey  et  al. ,  1986).  For  a  monotonically  increasing  stress  state,  the  current  model  is: 

e,  =  In (pjp)  =  (A  +  BaocJ(l-exp(-Caoel))  +  D(Toct/aocs)  [11 

where  e,  —  natural  volumetric  strain,  ln(v/Vi) 

p,  v  =  bulk  density  and  specific  volume  at  stress  state  and 
Pi,  Vj  =  initial  uncompacted  virgin  bulk  density  and  specific  volume 
t«a  —  applied  octahedral  normal  and  shear  stresses 
A,  B,  C  =  compactibility  coefficients 

D  =  coefficient  for  the  component  of  natural  volumetric  strain  due  to  applied  octahedral 
shearing  stress. 

The  two  idealized  boundary  conditions  of  (1)  zero  strain  at  a  zero  stress  state  and  (2)  linearly  asymptotic 
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at  high  hydrostatic  stress  states,  proposed  by  Bailey  et  al.  (1984,  1986),  are  maintained  in  this  model. 
This  model  (equation  1)  has  an  upper  bound  of  octahedral  shear  stress  at  which  plastic  flow  (strain  at 
constant  volume)  is  initiated  described  by: 

where  and  are  stress  values  at  maximum  density  and  K  is  a  coefficient  representing 

when  yield  is  initiated  by  plastic  flow.  This  follows  the  Drucker-Prager  failure  criteria  with  an  intercept 
of  zero  (Desai  and  Siriwardane  1984). 

Fig.  1  presents  bulk  density  data  from  two  Norfolk  sandy  loam  samples  displayed  as  a  function 
of  the  octahedral  stress  ratio,  The  two  straight  lines  representing  the  two  different  octahedral 

normal  stress  levels  have  the  same  slope.  The  evaluation  of  straight  line  slopes  at  all  levels  of  octahedral 
normal  stress  indicated  that  die  slope  was  independent  of  level  of  octahedral  normal  stress. 


Octahedral  stress  ratio 

Fig.  1.  Deviatoric  loading  portion  of  two  Fig.  2.  Typical  natural  shearing  strain  ratio 

Norfolk  sandy  loam  tests  as  a  function  of  the  data  for  Hiwassee  clay, 

octahedral  stress  ratio. 

A  shearing  strain  model  for  soil  that  includes  soil  behavior  under  compressive  normal  and  shear 
stresses  great  enough  to  attain  maximum  compaction  was  developed  (Johnson  and  Bailey,  1990). 
Representative  data  for  the  ratio  of  maximum  natural  shear  strain  to  the  volumetric  strain  occurring  after 
application  of  shear  stress  versus  the  ratio  of  maximum  shear  stress  to  major  principal  stress  are  shown 
in  Fig.  2.  The  maximum  natural  shearing  strain,  was  defined  as  die  difference  between  die  major 
and  minor  principal  natural  strains  according  to  Ludwik  in  1909  as  reported  by  Hoffman  and  Sachs 
(1953).  The  data  in  Fig.  2  are  for  six  stress  loading  paths  and  include  data  for  three  stress  loading  paths 

~  Clf  Cj,  C,)  from  Grisso  et  al.  (1987).  The  volumetric  strain  occurring  after  application  of 
shear  stress  in  the  ordinate  term  appears  to  account  for  much  of  die  stress  loading  path  or  stress  history 
effect.  These  data  (Fig.  2)  suggest  that  one  form  of  the  relationship  is: 

r«„/a,  =  K’  (1  -  0  exp(-h  Tw/e,  J)  [3] 

This  relationship  form  with  alternative  "variables"  of  and  which  appears  to  be  more 
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compatible  with  variables  in  equations  1  and  2,  is  being  investigated.  Four  soils,  Decatur  clay  loam 
(Rodic  Paleudults),  Hiwassee  sandy  loam  (Typic  Rhodulults),  Hiwassee  clay  (Typic  Rhodulults,  formerly 
classified  as  Lloyd  clay),  and  Norfolk  sandy  loam  (Typic  Paleudults)  from  the  NSDL  soil  bins  were  used 
to  develop  these  models. 

At  first  glance  these  model  results  may  appear  quite  different  than  'critical  state”  concepts  and 
resulting  cam  clay  and  cap  models  as  reported  in  the  literature  (Roscoe  el  al.  (19S8)  and  Desai  and 
Sirwardane  (1984))  that  were  developed  for  saturated  soils  using  the  concept  of  effective  stress. 

However,  equations  1  and  2  geometrically  represent  a  three-dimensional  surface,  illustrated  in  Fig.  3, 
which  is  similar  to  the  "critical  state"  surface  described  by  Roscoe  et  al.  (1958).  Fig.  4  illustrates  a 
form  of  a  cap  model  which  is  bounded  by  equation  2  and  equation  1  with  the  natural  volumetric  strain 
being  held  constant  (plastic  flow  condition)  at  a  value  caused  by  a  hydrostatic  stress  of  500  kPa  and 
t<Jo9 „  <.  K  (equation  2). 


Fig.  3.  Surface  represented  by  models  Fig.4  Projection  of  model  boundaries  on  the 
(Equations  1  and  2)  for  Hiwassee  clay.  octahedral  stress  plane  at  a  constant  volumetric 

strain  for  Decatur  clay  loam. 


Future  Plans 

Elastic  rebound  properties  will  be  determined  from  data  currently  being  collected  and  analyzed 
for  repeated  loading  and  unloading  in  a  conventional  triaxial  cell  using  two  different  stress  paths 
(constant  and  constant  cell  pressure)  (Johnson  et  al. ,  1992).  These  additional  elastic  rebound 
properties  and  tensile  stress-strain-strength  characteristics  are  needed  to  fully  implement  an  elasto-plastic 
or  "critical  state"-cap  type  finite  element  model  of  the  soil  behavior. 

The  role  of  the  intermediate  principal  stress  will  be  investigated  using  apparatus  designed  and 
developed  by  Gibas  et  al.  (1993).  This  will  test  the  validity  of  the  models  developed  from  use  of 
conventional  triaxial  cell  (an  axi symmetric  stress  state)  for  true  three-dimensional  stress  states. 
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FINITE  ELEMENT  MODELING  OF  WHEEL  PERFORMANCE 

AND 

SOIL  REACTION  AND  DEFORMATION 


Clarence  E.  Johnson,  Winfred  A.  Foster,  Jr. ,  Sally  Shoop  and  Randy  L.  Raper1 

Introduction 

Operation  of  wheeled  vehicles  over  the  land,  off  of  developed  roadways,  is  vital  to 
the  national  security  and  economy  of  the  United  States.  For  example,  part  of  our  national 
defense  relies  on  transport  of  supplies  and  manpower  provided  by  wheeled  vehicles  in 
regions  where  there  are  minimal  or  no  roadways.  Also,  our  mechanized  agricultural  and 
forest  industries  depend  on  cost  effective  and  efficient  wheeled  vehicle  and  equipment 
operation  on  soil  without  roadways  to  produce  food  and  fiber  for  our  national  economy. 

Thus,  it  is  important  to  develop  technology  to  predict  tractive  performance  of  wheeled 
vehicles  and  the  reaction  of  soil  to  wheeled  vehicle  traffic  under  a  variety  of  wheel 
configurations  and  soil  conditions.  This  technology  will  aid  the  design,  development  and 
utilization  of  future  wheeled  vehicles  with  improved  efficiency,  effectiveness  and/or  economy 
without  adverse  environmental  impact.  This  research  project  was  initiated  within  the  past 
year  and  has  the  following  progressive  objectives  as  follows: 

1.  Develop  a  plane  strain  finite  element  model  for  the  analysis  of  a  rigid  wheel  rolling 
on  the  edge  of  a  semi-infinite  linearly  elastic  plane. 

2.  Develop  a  plane  strain  finite  element  model  for  the  analysis  of  a  rigid  wheel  rolling 
on  the  edge  of  a  semi-infinite  elasto-plastic  plane. 

3.  Develop  a  plane  strain  finite  element  model  for  the  analysis  of  a  rigid  wheel  rolling 
on  tiie  edge  of  a  semi-infinite  plane  of  soil. 

4.  Expand  the  objectives  1,  2  and  3  to  a  full  three-dimensional  analysis. 

5.  Expand  the  objectives  1,  2,  3  and  4  to  include  a  non-rigid  wheel  where  part  or  all  of 
the  wheel  could  be  considered  "elastic". 


‘The authors  are:  Clarence  E.  Johnson,  Professor,  Agricultural  Engineering  Dept.,  Alabama 
Agricultural  Experiment  Station,  Winfred  A.  Foster,  Jr.,  Aerospace  Engineering  Dept.,  Aubum 
University,  AL.,  Sally  Shoop,  Research  Civil  Engineer,  U.S.  Army  Corp  of  Engineers,  Cold 
Regions  Research  and  Engineering  Laboratory,  Hanover,  NH  and  Randy  L.  Raper,  Agricultural 
Engineer,  National  Soil  Dynamics  Laboratory,  USDA-ARS,  Aubum,  AL. 
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Procedure 


All  models  defined  in  the  objectives  will  be  analyzed  uring  ABAQUS  and  some  will 
also  be  analyzed  using  NASTRAN  (finite  element  software)  for  comparison  purposes.  Data 
collected  by  Block  (1991)  at  the  National  Soil  Dynamics  Laboratory  (NSDL),  Aubum,  AL  in 
two  soils  with  an  instrumented  powered  rigid  wheel  will  be  utilized  to  validate  and  calibrate 
the  models.  These  data  include  data  from  five  pressure  cells  spaced  across  the  wi'ith  of  the 
wheel,  data  from  force  transducers  measuring  normal  force  and  tangential  force  on  the  face 
of  the  wheel  in  contact  with  the  soil,  and  stress  state  data  within  the  soil  at  two  depths 
beneath  the  path  of  the  wheel.  Data  from  various  kinds  of  triaxial  tests  for  development  of 
the  NSDL-AU  constitutive  model  parameters  (Bailey,  Johnson  and  Schafer  (1986)  and  Bailey 
and  Johnson  (1989))  are  also  available  for  these  two  soils. 

Constitutive  models  of  soil  behavior  for  objectives  3  and  4  similar  to  the  modified 
Cam  Clay,  Critical  State  and  the  NSDL-AU  constitutive  models  may  be  utilized. 

Current  Status 

As  a  starting  point,  we  decided  to  simulate  some  circular  plate  (approx.  18-in  dia)  and 
spherical  body  sinkage  data  that  Raper  (1987)  had  collected  in  the  NSDL  soil  bins.  Data  of 
force  vs  sinkage  and  soil  stress  state  at  four  locations  are  available.  A  linear  elastic 
axisymmetric  model  in  both  NASTRAN  and  ABAQUS  for  the  circular  plate  situation  using 
approximately  the  same  grid  size,  etc.  that  Raper  had  used  was  developed.  This  would  allow 
us  to  easily  make  comparisons  with  results  from  his  program  also.  Deformation  loading  with 
a  circular  plate,  without  gravitational  loading,  gave  "same"  results  from  all  three  finite 
element  programs. 

Gravitational  loading  (a  stress  boundary  condition)  combined  with  deformation  loading 
(a  "geometry"  boundary  condition)  presented  problems  in  both  NASTRAN  and  ABAQUS. 

So  next,  we  modeled  the  soil  being  loaded  with  a  "massless"  steel  circular  plate  on  the  soil 
surface  with  vertical  stresses  acting  near  the  center  of  the  plate.  This  approach  alleviated  the 
"combined"  boundary  problem  yet  allows  the  soil  to  experience  a  deformation  like  loading 
since  the  plate  is  very  rigid  compared  to  the  soil. 

We  found  that  the  combination  of  non-linear  elasticity  and  axisymmetry  requires  full 
3D  model  elements.  So  a  3D  grid  for  an  axisymmetric  section  was  developed  that  could  be 
used  in  both  ABAQUS  and  NASTRAN.  Currently,  we’re  using  a  non-linear  bulk  modulus 
(tangent)  data  array  we  developed  from  hydrostatic  soil  compaction  data  in  a  triaxial  cell  for 
the  soil  used  in  the  "plate  sinkage  tests"  and  a  constant  Poisson’s  ratio  of  0.35.  Results 
from  the  linear  elastic,  nonlinear  elastic,  and  elasto-plastic  behavior  models  have  similar 
mean  normal  stress  distributions  in  the  soil  which  compare  favorably  with  Raper ’s  data  in  a 
Norfolk  sandy  loam  soil.  Displacements  within  the  soil  are  under  predicted  by  all  three 
behavior  models. 

Future  Plans 
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Our  next  step  is  to  develop  interface  or  contract  elements  to  work  with  a  curved 
surface.  These  are  the  elements  between  a  spherical  surface  and  the  soil  surface  in  an 
axisymmetric  model  or  between  a  cylindrical  surface  (rigid  "wheel")  and  the  soil  surface  in 
plane  strain  and  plane  stress  models.  This  should  help  us  when  using  either  deformation  or 
"stress"  loading  with  a  "rigid"  curved  surface  at  the  soil  surface. 

References 

Bailey,  A.C.,  C.E.  Johnson,  and  R.L.  Schafer  (1986)  "A  model  for  agricultural  soil 
compaction " ,  J.  agric.  Engng  Res.  33:257-262. 

Bailey,  A.C.  and  C.E.  Johnson  (1989)  "A  soil  compaction  model  for  cylindrical  stress 
states",  TRANSACTIONS  of  the  ASAE  32(3):  822-825. 

Block,  W.A.  (1991)  Analysis  of  soil  Stress  under  rigid  wheel  loading  Geologic 
Materials,  Unpublished  Ph.D.  Dissertation.  Auburn  University,  Auburn,  AL. 

Raper,  R.L.  (1987)  Prediction  of  soil  compaction  using  the  finite  element  method. 
Unpublished  Ph.D.  Dissertation.  Iowa  State  University,  Ames,  IA. 


3 


GENERALIZED  JANOSES  SHEAR  STRESS-SLIPPAGE  RELATION 

Hidenori  Murakami1  and  Tatsunori  Katahira1 


ABSTRACT 

The  authors  propose  a  shear  stress-slippage  relation  for  plane  slippage  in  arbitrary 
directions  to  furnish  a  soil-track  interaction  relation  for  off-road  mobility  analyses  of  tracked 
vehicles.  The  relationship  is  obtained  based  upon  a  suite  of  plate  shear  tests  conducted  in  a  soil  bin 
of  loose,  dry  sand,  and  generalizes  Janosi's  shear  stress-slippage  relation  proposed  for  slippage  in 
longitudinal  as  well  as  lateral  directions. 

DESCRIPTION  OF  EXPERIMENTS 

In  order  to  develop  a  model  for  shear  stress-slippage  relation  for  tracks,  plate  shear  tests 
were  conducted  in  a  soil  bin  filled  with  dry,  loose  sand  (Fig.  la).  Two  types  of  plates  --  with  and 
without  grousers  --  of  width  8  cm  and  length  42  cm  were  tested.  The  grouser  pitch  is  2  cm  and  the 
height  is  0.7  cm  (Fig.  lb).  A  special  load  cell  was  employed  to  measure  two  shear  force 
components  under  prescribed  slippage,  slip  velocity,  and  normal  force.  The  slip  angle,  0,  is 
measured  from  the  longitudinal  (xr)  axis  of  the  track  in  the  clockwise  direction  in  the  plan  view. 
For  a  prescribed  set  of  slip  direction,  slip  velocity,  and  normal  force,  the  shear  force  components 
in  the  direction  of  slippage,  Qs,  and  in  the  transverse  direction,  Qt,  were  measured  along  with  the 
sinkage. 


. .  .  vertical  load 
driving  motor 


sinkage 

sensor 


slippage  sensor 


Fig.  la  The  shear  test  apparatus 
EXPERIMENTAL  DATA  AND  ANALYSES 


Experimental  Results 

The  shear  force  components  and  sinkage  were  measured  for  monotonic  slippage  in  the  slip 
directions,  0  =  0°  (longitudinal  slippage)  to  90°  (lateral  slippage)  for  every  15°at  constant  slip 
velocity,  0.3  cm/s.  For  a  plate  with  grousers,  the  shear  force  components,  Qj  versus  slippage 
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and  Qt  versus  slippage,  for  the  same  slip  directions  are  illustrated,  respectively,  in  Figs.  2a  and  2b. 
Similar  results  were  obtained  for  a  plate  without  grousers. 


slippage  [us]  (mm) 


Fig.  2a  Shear  force  component  0$  versus  slippage  for  a  flat  plate  with  grousers 
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Fig.  2b  Shear  force  component  Qt  versus  slippage  for  a  flat  plate  with  grousers 


For  longitudinal  and  lateral  directions  (0=0*  and  0=90’  respectively)  the  shear  force 
component,  Qt,  vanishes  as  shown  in  Fig.  2b.  The  monotonic  loading  curves  in  the  longitudinal 
and  lateral  directions  were  described  by  Janosi  and  Hanamoto  (1961)  as 


Qeff  =  Qy(Q3){  1-exp  {-ptUefflH  (1) 

where  is  the  shear  force  in  the  slip  direction,  Qy  is  the  critical  shear  force,  [u^]  is  the  slippage, 
and  p  is  a  constant  The  critical  shear  force  Qy  increases  with  increasing  normal  force  Qj. 

For  other  slip  directions  Qt  changes  with  slip  directions.  The  comparisons  of  Figs.  2  with 
the  results  without  grousers  have  revealed  that  indeed,  the  plate  with  grousers  exhibits  orthotropic 


dependency  of  shear  force-slippage  relations.  The  objective  of  this  paper  is  to  establish  appropriate 
definitions  of  and  [ueff]  so  that  all  curves  in  Figs.  2  can  be  deduced  from  a  single  master 
curve. 

Slip  Surface  in  die  Interaction  Shear  Force  Plane 

In  order  to  find  the  shape  of  the  critical  slip  surface  in  the  interaction  shear  force  plane,  the 
shear  forces  Q,  and  Qt  are  transformed  into  the  track  axial  and  lateral  components,  and  Q2, 
according  to  the  coordinate  transformation  between  the  xlt  x2  coordinate  system  attached  to  the 
track  and  the  s,  t  coordinate  system  which  described  the  slippage  and  transverse  directions.  The 
slippage  [u]  is  also  decomposed  into  the  xt  and  x2  components,  [ut]  and  [uj. 


Qi(N)  Q,(N) 


Figs.  3  Loading  Surfaces  for  the  loading  plate  with  grousers  under  Q3=98N  and  196  N 

Figures  3  show  the  loading  surfaces,  at  plastic  slippage  7  mm,  17  mm,  and  27  mm,  in  the 
shear  force  plane  for  the  plate  with  grousers,  for  the  normal  force  Q3  =  98N  and  196  N.  The 
experimental  data  for  both  Qj  =  98N  and  196  N  show  that  loading  surfaces  in  the  Qt  and  Q2  plane 
can  be  described  by  ellipses.  The  arrows  of  plastic  slip  velocity  plotted  on  the  loading  surface 
show  that  plastic  slip  velocity  is  normal  to  the  loading  surface. 

Generalized  Janosi's  Shear  Force-Slippage  Relation 

For  an  arbitrary  plane  slippage,  the  effective  shear  force  Qeff  is  introduced: 

Qrfr  =  V(Qi)J  +  (KQ2)2  .  (2) 

where  k  represents  the  ratio  between  the  longitudinal  critical  shear  force  and  the  lateral  critical  shear 
force.  Equation  (2)  describes  an  elliptical  loading  surface  in  the  interaction  shear  force  plane. 

In  order  to  account  for  irreversible  slippage  observed  after  unloading,  the  slip  velocity  is 
decomposed  into  elastic  and  plastic  parts  denoted  by  fVj]®1  and  [VJP1.  From  the  normality  of  the 
slip  velocity  to  each  loading  surface  the  following  effective  slip  velocity  is  employed: 


Figures  4  show  the  data  in  Figs.  2  expressed  with  respect  to  the  above  effective  quantities. 

The  effective  shear  stress  is  defined  as  xeff=  QefJ/A  where  A  is  the  contact  area  of  the  loading  plate. 
The  results  show  that  all  the  curves  in  Figs.  2  collapse  nicely  with  the  effective  shear  force  and 
slippage,  and  the  collapsed  loading  curves  are  described  by  Janosi's  monotonic  loading  curve  (1). 
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slippage  [Ujff]  (mm) 


Figs.  4  The  effective  shear  force-slippage  relation  for  the  loading  plate 
with  grousers  under  Q3=98N  and  196  N 
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Introduction 


At  KRC  our  main  area  of  mobility  research  and  development  has  been  in  off-road  terrains  in 
cold  climates  which  includes  shallow  snow,  deep  snow,  thawing  soils,  compacted  snow  and 
ice.  We  have  done  some  modeling,  field  validation  and  development  of  devices  for 
characterizing  the  different  terrain  materials.  The  following  is  a  brief  summary  of  the  areas 
that  we  have  been  involved  with  in  the  past. 

Modeling  -  KRC  has  run  the  NRMM  for  several  clients  over  the  past  several  years  and  have 
worked  with  CRREL  and  WES  in  the  development  and  validation  of  the  Shallow  Snow 
Model  for  the  NRMM..  We  have  developed  a  simple  vehicle  countermobility  model  in  snow. 
KRC  also  has  a  group  that  has  developed  a  sophisticated  thermal  signature  model  of  the 
terrain  and  vehicles  which  has  been  widely  distributed  among  Army  and  Air  Force  users. 
Overall,  KRC  has  had  a  minimal  amount  of  experience  in  developing  new  mobility  models. 

Field  Validation  and  Evaluation  -  KRC  has  a  tremendous  amount  of  experience  in  evaluating 
vehicles  for  mobility  performance  in  winter  terrains  (and  validating  model  predictions).  We 
have  tested  many  military  tracked  vehicles  ranging  from  a  snowmobile  up  to  the  M1A1  tank. 
KRC  has  also  evaluated  several  wheeled  vehicles  ranging  in  size  from  small  jeeps  up  to  5  ton 
trucks.  We  jointly  participated  in  the  "Wheels  vs  Tracks"  project  with  CRREL  and  WES. 
KRC  has  been  involved  in  unique  projects  such  as  a  comparison  between  a  standard  wheeled 
MK48  and  a  prototype  version  utilizing  the  Caterpillar  Mobil  Trac  System.  We  have  tested 
vehicles  with  anti-lock  brakes,  various  traction  systems  and  CTIS.  Types  of  tests  conducted 
in  the  past  have  included  drawbar  pull,  motion  resistance,  mobility  evaluations,  slope 
climbing,  side  slope  evaluation,  acceleration,  braking  and  handling.  Although  the  majority  of 
our  testing  experience  has  been  in  winter  terrains,  especially  snow  and  ice,  the  multi  season/ 
terrain  capability  of  several  vehicles  has  been  evaluated.  The  wide  range  of  vehicles,  terrains, 
and  tests  has  resulted  in  KRC  developing  expertise  in  instrumentation,  data  acquisition,  data 
processing  and  validation.  Current  equipment  includes  a  wide  range  of  sensors  and 
transducers,  portable  computerized  data  acquisition  systems  and  telemetry  systems.  Finally, 
our  work  in  the  winter  mobility  area  has  resulted  in  the  design,  development  and  testing  of 
ice  cleats  for  the  Marine  Corps  LTV-P7  and  the  M1A1  Abrams  tank. 


Terrain  Characterization  -  KRC  has  performed  various  soil  property  measurements  including 
use  of  the  cone  penetrometer  but  mobility  in  soils  has  not  been  our  main  area  of  interest. 

Snow  is  another  story.  KRC  was  one  of  the  first  to  build  a  bevameter  including  both  the 
shear  and  compaction  device.  We  have  used  several  snow  density  kits  and  developed  one  of 
our  own.  We  have  also  used  the  Ramsonde,  Canadian  snow  hardness  gauges,  several  types  of 
free  water  content  methods  and  coefficient  of  friction  or  traction  devices  for  characterization 
of  snow.  KRC  has  designed  and  built  a  portable  bevameter  as  well.  Currently,  KRC  is 
testing  a  CRREL  developed  load  frame  device  for  characterizing  snow.  On  ice,  we  have  used 
a  friction  tester  for  measuring  the  coefficient  of  friction  of  the  ice.  KRC  has  also 
conceptualized  a  design  for  an  ice  strength  measurement  device  but  has  never  fully  designed 
or  fabricated  one.  Thus  far,  we  have  not  been  able  to  find  any  other  organization  interested 
in  measuring  the  strength  of  ice.  As  far  as  KRC  is  concerned,  when  we  develop  a  traction 
aid  for  a  tracked  vehicle  and  test  two  or  more  designs  on  two  or  more  different  days  we  need 
to  know  the  characteristics  of  the  ice.  Ice  tends  to  change  and  coefficient  of  friction  is  not 
useful  because  the  ice  cleats  dig  in  to  the  ice.  Penetration  and  shear  strength  are  the 
important  parameters  with  traction  devices  for  ice. 

Future  Areas  Of  Interest 


We  are  currently  working  with  vehicles  that  have  anti-lock  brakes,  various  traction  control 
systems,  central  tire  inflation  systems  and  independent  suspensions.  It  has  always  appeared 
that  people  have  considered  vehicle  dynamics  modeling  separate  from  vehicle  mobility 
modeling.  At  KRC,  we  would  like  to  bring  the  vehicle  dynamics  closer  to  mobility  modeling 
for  modeling  ABS  and  TC  systems  in  off-road  winter  terrains.  This  will  require  being  able  to 
input  deformable  terrains  into  a  dynamics  model.  A  better  tire  model  is  also  required.  Some 
tire  models  are  known  but  most  are  not  considered  ideal.  We  know  that  some  companies, 
i.e.,  tire  companies  and  some  of  the  ABS  and  auto  companies  have  better  tire  models  but  they 
don’t  want  to  release  them  because  they  have  put  all  their  resources  into  it  and  consider  them 
proprietary  information.  "Good"  tire  models  for  hard  pavement  may  be  insufficient  for  off¬ 
road  modeling.  Our  future  efforts  will  include  attempting  to  use  the  models  of  deformable 
terrain  in  the  dynamics  models. 


USING  THE  FINITE  ELEMENT  METHOD  TO  PREDICT  SOIL 
STRESSES  BENEATH  A  RIGID  WHEEL 

R.L.  Raper1,  C.E.  Johnson2,  A  C.  Bailey1,  and  E  C.  Burt1 


Introduction 

The  objective  of  this  experiment  was  to  investigate  the  ability  of  the  finite  element 
method  to  predict  soil  stresses  beneath  a  rigid  wheel  in  two  soils  and  in  two  soil  conditions 
An  experiment  was  conducted  in  the  soil  bins  at  the  National  Soil  Dynamics  Laboratory 
(NSDL)  during  which  the  transducers  were  used  to  measure  soil  stress  beneath  the  rigid 
wheel.  Two  different  constitutive  relationships  for  soil  were  compared  to  determine  which 
modeled  actual  soil  behavior  the  closest.  Soil  stress  measurements  were  then  compared  to 
results  predicted  with  the  finite  element  method. 


Procedure 

The  plane  strain  assumption  was  used  to  model  the  rigid  wheel.  The  rigid  wheel  must 
be  visualized  as  infinitely  wide  to  understand  this  assumption.  This  of  course  is  not  true,  but 
the  stresses  beneath  the  center  of  the  rigid  wheel  should  not  differ  greatly  from  those  beneath 
an  infinitely  wide  cylinder. 

Modeling  the  soil  matrix  is  the  most  difficult  problem  being  confronted  by  soil 
compaction  researchers.  This  non-homogeneous,  non-linear,  elastic-plastic,  particulate 
medium  makes  exact  solutions  impossible.  Assumptions  must  be  made  about  its  previous 
history,  the  existence  of  clods,  and  the  location  of  hard  pans.  Neglecting  these  problems  and 
treating  the  soil  as  a  homogeneous  medium  has  allowed  some  limited  successes  in  soil 
compaction  modeling.  A  model  has  been  developed  that  relates  the  volumetric  strain  to  the 
applied  hydrostatic  stress  (Bailey  et  al,  1984). 

ev  =  (A+Bo^*(  l-ei~Co*])  O) 

where  e¥  =  volumetric  strain,  (change  in  volume  /  original  volume) 
ahyd  =  hydrostatic  stress,  kPa 

A,  B,  and  C  =  compatibility  coefficients  established  by  fitting  data  to  equation. 
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The  hydrostatic  soil  compaction  model  assumes  that  all  stresses  surrounding  a  soil 
particle  exerted  equal  forces  on  this  particle.  Of  course  this  assumption  does  not  allow  for 
the  development  of  unequal  directional  stresses  which  create  shear.  These  shear  stresses  have 
been  shown  to  increase  the  magnitude  of  soil  compaction.  A  model  has  been  developed  that 
improves  on  the  hydrostatic  soil  compaction  model  to  include  the  effects  of  shear  stress 
(Bailey  and  Johnson,  1990). 


e,  *  ) 

o _ , 


(2) 


where  e,  =  natural  volumetric  strain,  In  (change  in  volume  /  original  volume) 

Oojt  =  octahedral  normal  stress,  kPa 
=  octahedral  shearing  stress,  kPa 
D  =  another  compactibility  coefficient 

This  model  is  similar  to  Equation  1  except  for  the  addition  of  the  shearing  stress  component 
and  the  use  of  the  natural  strain  definition.  Limitations  were  placed  on  the  shearing  stress 
component  in  the  above  model  to  indicate  maximum  density  at  plastic  flow.  The  restricting 
relationship  is: 

X  =  K*~  (3) 

where  to,.  =  ultimate  shearing  stress  at  maximum  density 
K  =  coefficient  representing  soil  plastic  flow  yield 

Although  the  soil  constitutive  equations  have  been  given,  much  more  must  be 
accomplished  before  they  are  useable  in  a  finite  element  model.  These  equations  must  be 
used  to  predict  the  linear-elastic  parameters.  Young's  Modulus  (E)  and  Poisson's  ratio  (v). 

Each  of  these  parameters  varies  with  different  stress  and  strain  levels  in  the  soil  (Duncan  and 
Chang,  1970).  For  further  information  regarding  the  prediction  of  these  parameters  using 
triaxial  data  see  Raper  et  al.  (1992) 

The  mesh  designed  to  model  the  soil  beneath  the  rigid  wheel  is  shown  in  Figure  1 
along  with  a  deformed  grid  overlaid.  The  exact  shape  of  the  rigid  wheel  was  maintained  as  it 
came  in  contact  with  the  soil  until  it  applied  the  maximum  load.  Twenty  load  steps  were 
found  to  be  adequate  to  allow  the  model  to  incorporate  the  non-linear  behavior  of  soil. 

An  experiment  was  performed  using  the  rigid  wheel  in  the  soil  bins  at  the  NSDL. 

This  wheel  (30.5  cm  wide  and  137.2  cm  in  diameter)  was  used  in  two  different  soil  types,  a 
Norfolk  sandy  loam  soil  and  a  Decatur  clay  loam  soil.  The  soils  were  prepared  in  a 
uniformly  loose  state  and  also  with  a  hard  pan.  The  experiment  as  reported  here  involved 
dynamic  loads  on  the  rigid  wheel  of  5.8  and  11.6  kN.  Four  replications  were  performed  of 
each  treatment. 

Stress  State  Transducers  (SSTs)  (Nichols  et  al,  1987)  were  buried  beneath  the  center 
of  the  rigid  wheel  at  a  30  cm  depth  in  the  loose  soil  and  on  the  hard  pan  in  the  other 
treatment.  The  final  depths  of  the  SSTs  were  used  to  establish  the  depth  that  the  finite 
element  results  would  be  analyzed. 


Figure  1.  Original  finite  element  mesh  and  final  displaced  mesh  showing  location  of 
impeding  layer  in  Decatur  clay  loam  soil  when  loosely  tilled.  Only  one  half  of  the  soil 
was  modeled  beneath  die  rigid  wheel  because  of  symmetry. 


Displacements  of  the  soil  surface  and  the  transducers  were  measured  at  the  conclusion 
of  each  experimental  run.  The  surface  displacement  depths  were  used  to  load  die  finite 
element  model  in  the  vertical  direction.  An  average  depth  of  surface  displacement  was 
obtained  for  each  treatment  and  a  95%  confidence  interval  established.  Three  finite  element 
models  were  then  run;  one  at  the  mean  value,  and  one  each  at  die  upper  and  lower  95% 
confidence  intervals. 

Penetrometer  measurements  were  made  in  die  soil  bins  to  determine  die  depth  of  any 
impeding  layers.  In  the  hard  pan  treatments,  this  layer  was  found  at  a  depth  of  approximately 
36  cm  in  both  die  Norfolk  sandy  loam  soil  and  die  Decatur  clay  loam  soil.  When  die  loose 
soil  treatment  was  used,  a  root  impeding  layer  was  found  at  a  depth  of  48  cm  in  die  Norfolk 
soil  and  at  54  cm  in  the  Decatur  soil.  These  depths  were  used  to  fix  the  nodes  of  the  finite 
element  mesh  to  prevent  soil  movement  past  this  depth. 


Results  and  Discussion 

The  peak  octahedral  normal  stress  and  peak  major  principal  stress  were  investigated  to 
determine  if  they  fit  within  95%  confidence  intervals  of  stress  determined  from  die  SSTs. 
When  examining  die  octahedral  normal  stresses  (Figure  2),  the  hydrostatic  model  seems  to  be 
the  better  model  if  only  the  lower  loads  are  considered.  This  result  may  be  reasonable 
because  of  the  lack  of  shear  stress  that  is  developed  at  the  low  load  levels,  which  the 
hydrostatic  model  does  not  account  for.  When  considering  die  high  load  treatment,  however, 
neither  model  was  able  to  predict  the  stresses  with  much  certainty,  especially  in  the  Norfolk 
soil. 

Major  principal  stress  was  predicted  with  slighdy  more  accuracy  (Figure  2).  Again  at 
the  low  load  levels,  die  hydrostatic  model  fit  across  both  soil  types  and  both  soil  conditions. 
At  the  high  loads,  the  shear  stress  model  fit  all  of  the  measured  data  except  the  high  load 
treatment  in  the  Norfolk  soil.  The  shear  stress  model  managed  to  fit  all  of  die  data  in  the 
Decatur  soil,  both  low  and  high  loads. 


Figure  2.  Octahedral  normal  stress,  major  principal  stress,  o,,  and  their  95%  confidence 
intervals  measured  with  the  SSPs  plotted  against  die  finite  element  results. 

Continued  development  of  the  finite  element  model  is  warranted  to  allow  better 
predictions  of  soil  stress  to  be  accomplished.  A  true  three-dimensional  finite  element  model 
could  enable  better  predictions  to  be  made.  Excessive  stress  predictions  of  die  shear  stress 
model  could  be  the  result  of  the  plane  strain  assumption.  These  large  stresses  could  be  due  to 
the  effect  of  confining  stresses  that  develop  from  the  rigid  wheel  being  modeled  as  an 
infinitely  long  roller. 
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A  Contact  Mechanics  Approach  to  the  Modeling  of 
Dynamic  Soil-Vehicle  Interaction 

Dr.  Antoinette  TordesiUas1 


Introduction 


A  new  approach  to  the  modeling  of  soil-vehicle  interaction  is  introduced.  The 
interaction  at  the  interface  between  vehicle  and  soil  is  formulated  as  a  dynamic  contact 
problem  and  is  solved  using  the  principles  and  methodologies  of  the  theory  of  Contact 
Mechanics.  This  approach  has  two  important  advantages  over  existing  analytical  soil- 
vehicle  interaction  models.  First,  a  contact  mechanics  formulation  obviates  the  need  to 
know  a  priori  the  stresses  or  deformations  at  the  soil-vehicle  interface.  Instead,  these 
interfacial  properties  are  determined  using  directly  and  precisely  measurable  quantities. 
Second,  this  formulation  avails  the  analysis  of  the  interfacial  phenomena  to  the 
comprehensive  theories  of  contact  mechanics  and  tribology,  with  their  proven  economic 
and  reliable  techniques  for  establishing  detailed  information  on  contact  properties.  The 
scope  of  the  field  of  contact  mechanics  is  extensive.  Material  models  which  have  been 
commonly  studied  span  the  range  from  elastic,  viscoelastic,  elastoplastic,  to  viscoplastic. 
In  such  analyses,  various  geometric  properties  and  contact  configurations  of  the  bodies 
have  been  considered,  in  conjunction  with  both  non-classical  and  classical  Coulomb 
friction  laws.  We  arc  conducting  a  preliminary  study  on  the  soil-tire  interaction  system.  In 
accordance  with  the  studies  of  Pi  (1988),  the  basic  soil  behavior  under  dynamic  vehicle 
passage  is  considered  to  be  viscoelastic  and  is  represented  by  a  Maxwell-Kelvin  3- 
parameter  model.  A  new  model  for  the  tire  is  introduced  which  consist  of  a  three- 
dimensional  circular  elastic  cylinder,  and  is  based  on  recently  derived  stress-displacement 
constitutive  relations  unique  to  the  cylindrical  geometry.  This  is  a  significant  improvement 
to  the  previously  adopted  Hertz  theory  in  which  the  cylinder  is  idealized  and  assumed  to 
deform  as  an  elastic  half-space.  Thus,  this  new  tire  model  incorporates  the  pertinent  tire 
curvature  and  edge  effects  into  the  overall  soil-tire  interaction  model. 

The  Theory  of  Two-Body  Contact  Mechanics 


The  theory  of  contact  mechanics  concerns  itself  entirely  with  the  local  interaction 
phenomena  arising  at  the  interface  between  two  bodies  which  are  brought  into  contact. 
This  emphasis  on  the  contact  interface  rests  on  the  premise  that  the  conditions  therein 
determine  the  internal  states  of  each  body.  Specifically,  if  the  fundamental  properties 
consisting  of  the  contact  stresses  and  deformations,  as  well  as  the  size  and  shape  of  the 
contact  area  are  known,  then  the  stresses  and  displacements  at  any  point  inside  each  body 
can,  in  principle,  be  found.  The  general  laws  of  contact  mechanics  arc  summarized  below 
and  are  illustrated  in  Figure  1: 

1)  Compatibility  condition:  no  interpenetration  exist  between  the  contacting  surfaces,  i.e. 


uzi  +ul2+h(x,y)-S 


=  0,  inside  £2, 
>  0,  outside  £1 


(1) 


where  S  is  the  relative  approach  of  the  bodies,  £2  is  the  contact  area,  h(x,y)  is  the  initial 
separation  of  surface  points,  and  uri(body  i=l,2)  denote  the  displacements. 

2)  No  tensile  tractions  are  allowed  on  contacting  boundaries.  The  normal  stress  p  is 
compressive  and  vanishes  outside  the  contact  area  £2, 

p  >  0,  inside  £2;  p-  0,  outside  £2.  (2) 
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Figure  1.  The  contact  of  two  deformable  bodies:  (a)  unloaded,  (b)  loaded. 

The  tangential  stresses  q  are  related  to  the  normal  stresses  by  the  appropriate  friction  law 
depending  on  the  nature  of  the  contact  surfaces.  If  Coulomb's  law  is  assumed,  then  two  zones  of 
adhesion  and  slip  within  the  contact  area  result,  viz.  and  £2S,  respectively,  and 

\q\  <  pip,  adhesion  zone  ClA;  q  =  ±pip,  slip  zone  £2S;  ClA  +  £ls  =  £1  (3) 

Solution  techniques  for  various  classes  of  contact  problems  have  been  developed  which 
incorporate  either  finite  element  or  boundary  element  methods.  The  finite  element,  in  conjunction 
with  a  variational  formulation  of  the  contact  system,  has  proved  to  be  the  most  powerful 
methodology  for  solving  the  more  difficult  classes  of  contact  problems  such  as  three-dimensional 
frictional  contact  The  premise  of  the  variational  theory  for  solid  body  contact  is  that  the  true 
contact  area  and  the  stresses  which  act  therein  are  those  which  minimize  an  appropriate  energy 
function.  On  this  basis,  a  set  of  variational  inequalities  which  constitute  a  minimization  problem  of 
functionals  can  be  derived.  Specifically,  the  contact  problem  can  be  reduced  to  a  single  relation  in 

which  the  total  complementary  energy  U*  =  f(u6,p,q,Q)  is  minimized  subject  to  the  conditions  of 
equations  (2)  and  (3).  In  a  numerical  implementation  the  resulting  minimization  problem  are 
discretized  to  yield  a  system  of  equations  which  can  be  solved  using  any  of  the  existing 
mathematical  programming  techniques  of  optimization.  Clearly,  the  distinct  advantage  of  this  type 
of  formulation  is  that  all  boundary  conditions,  including  the  contact  conditions,  are  incorporated 
into  a  single  variational  inequality  and  available  minimization  routines  can  be  used  to  solve  the 
problem.  Expected  results  from  a  contact  mechanics  model  include:  size  and  shape  of  the  contact 
area,  complete  surface  and  subsurface  stresses  and  displacements,  boundary  bordering  regions  of 
slip  and  adhesion  at  the  contact  interface,  and  the  relative  approach  of  the  bodies  (i.e.  sinkage). 

Constitutive  Stress-Displacement  Relations 

Of  prime  importance  in  a  contact  mechanics  formulation  are  the  body  displacements  uz 
which  result  from  the  surface  tractions  at  the  contact  interface,  as  called  for  in  the  compatibil.  v 
relation  of  equation  (1).  It  is  therefore  important  that  the  stress-displacement  relation  accurately 
represent  the  body's  response  to  surface  tractions,  from  the  point  of  view  of  both  its  material  and 
geometric  properties. 


SOIL  Mass.  The  stress-strain  behavior  of  soil  is  simulated  by  a  standard  three-parameter 
model  as  in  Pi  (1988),  consisting  of  a  Maxwell  spring  and  a  Kelvin  element  in  series  as  illustrated 
in  Figure  2.  It  has  been  shown  that  an  equivalent  of  the  Boussinesq-Cemiti  linear  elastic  relation 
can  be  derived  for  the  viscoelastic  half-space  using  the  correspondence  principle  of  elasticity 
(Kalker  1990).  This  facilitates  the  calculation  of  the  entire  elastic  field,  viz.  the  stresses  and 
deformations  both  on  the  surface  and  in  the  interior  of  a  viscoelastic  half-space. 


E2 


Vi 


Figure  2.  Viscoelastic  model  for  soil;  elastic  moduli  Ej  and  E2,  and  coefficient  of  viscosity  tj2- 

Tire.  For  the  tire  model,  we  introduce  a  three-dimensional  circular  elastic  cylinder  with 
diameter  and  width  equal  to  that  of  the  pneumatic  tire  to  be  analyzed.  The  constitutive  stress- 
displacement  relation  accounts  for  both  its  sectional  and  longitudinal  curvature  and  dimensions.  A 
recently  derived  solution  for  the  stress-displacement  relation  unique  to  a  circular  elastic  cylinder  in 
a  two-body  contact  system  is  adopted  (Tordesillas  and  England  1994).  This  should  provide  an 
improvement  on  models  which  are  based  on  the  Hertz  theory  in  which  the  curvature  of  the  bodies 
are  ignored.  Prior  to  our  work  in  Hill  and  Tordesillas  (1989,  1992),  various  simplifications  of  this 
type  were  adopted  to  the  study  of  cylindrical  contact  since  the  important  solution  for  a  point 
force(s)  acting  on  the  boundary  of  a  circular  elastic  cylinder  had  not  been  established.  Such  stress- 
displacement  solutions  relating  to  specific  point  force  systems,  as  shown  in  Figure  3(b),  yield  the 
constitutive  stress-displacement  relation  for  the  body  under  a  distributed  loading  as  depicted  in 
Figure  3(a).  The  procedure  is  based  on  the  classical  superposition  principle  of  linear  elasticity,  and 
involves  the  summation  of  the  point  force  solutions  over  the  contact  area  to  obtain  the 
corresponding  stress-displacement  relation  for  the  body  subject  to  a  distribution  of  forces  at  its 
boundary.  In  Hill  and  Tordesillas  (1989),  we  developed  a  novel  technique  to  derive  exact 
solutions  for  various  point  force  systems  relating  to  cylindrical  three-body  contact  problems.  The 
technique  is  based  on  complex  variable  theory  and  is  one  that  has  been  recently  used  to  establish 
the  basic  solution  corresponding  to  the  two-body  contact  problem.  The  empirical  input  parameters 
for  the  tire  model  are:  Diameter,  Width,  Poisson's  ratio,  vlin  ,  Young's  elastic  modulus,  Etire. 
Specifically,  the  Young's  elastic  modulus  for  the  tire  must  be  as  close  as  possible  to  that  for  the 
given  pneumatic  tire,  and  reflect  its  particular  inflation  pressure  and  carcass  stiffness  or  strength. 


(a)  (b) 


Figure  3.  (a)  Two-body  contact  involving  a  circular  elastic  cylinder,  (b)  Basic  point-force  system. 

The  phenomenon  of  stress  concentrations  arising  from  'edge  contact'  is  being  studied  for 
the  case  of  a  viscoelastic  media.  Contact  stress  concentrations  which  arise  from  any  discontinuities 


or  significant  changes  in  body  profiles  manifest  themselves  as  stress  singularities  when  studied 
within  the  framework  of  linear  elasticity.  A  thorough  evaluation  of  the  precise  structure  and  order 
of  these  singularities  was  carried  out  by  Comninou  (1976)  for  frictional  contact  of  bodies  of 
various  elastic  and  geometric  properties.  In  certain  cases,  the  appropriate  stress  concentrations  at 
the  edges  of  a  cylinder  can  be  incorporated  using  the  technique  employed  in  Tordesillas  and  Hill 
(1991).  This  technique  essentially  involves  scaling  the  discretized  stresses  at  each  element  i,  p(, 
by  the  appropriate  order  of  singularity.  This  procedure  has  been  successfully  applied  in  the  design 
analysis  of  roller  bearings  whose  edges  have  been  partially  rounded  at  t*e  ends  so  as  to  relieve 
stress  concentrations  (Ahmadi  et  al.  1983).  We  successfully  applied  it  to  the  analysis  of  cylindrical 
steel  and  rubber-covered  steel  contact  systems  (Tordesillas  1991). 
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Tire-Terrain  Modeling  for  Deformable  Terrain 

S.  Shoop,  CRREL 


Objective:  To  develop  a  numerical  model  simulating  the  interaction  between  a  tire  and 
highly  deformable  terrain  material.  Vehicle  movement  on  typical  cold  regions  surfaces 
may  result  in  large  deformation  from  both  compaction  and  mass  movement.  The  model 
will  be  of  assistance  in  predicting  vehicle/tire  performance,  designing  off-road  tires,  and 
estimating  terrain  damage  due  to  compaction  and  rutting. 

Approach:  Current  tire  design  technology  does  not  consider  the  interaction  of  the  tire 
and  various  deformable  materials  such  as  soil  or  snow,  focusing  primarily  on  tire/pavement 
interactions.  This  project  is  designed  to  integrate  research  being  conducted  in  two  very 
different  areas;  experimental  and  numerical  simulation  of  tractive  loading  on  deformable 
terrain,  and  numerical  models  of  tire  deformation. 

Two-dimensional  simulations  will  be  used  to  study  the  kinematics  of  the  problem.  Three- 
dimensional  simulations  will  concentrate  initially  on  the  terrain  deformation,  and  finally  the 
full  tire-terrain  system. 

Much  preliminary  work  has  been  completed  that  relates  directly  to  the  modeling  effort 
including  constitutive  models  of  cold  regions  terrain  (snow,  freezing  ground,  thawing 
ground),  detailed  measurement  of  terrain  deformations  including  melting,  and  numerous 
measurements  of  the  forces  generated  at  the  tire-surface  interface  for  off-road,  unsurfaced 
roads,  freezing  and  thawing  soils,  snow  and  ice. 

Applications: 

•  Parametric  study  of  effects  of  tire  configuration  and  snow/soil  properties  on  tire 
performance  and  terrain  deformation/damage. 

•  Effects  of  strong  soil  layering  (such  as  produced  by  freeze/thaw)  on  tire/vehicle 
performance 

•  Damage  prediction  on  unpaved  surfaces 

•  Off-road  and  all-season  tire  design 

Results  will  be  of  interest  to  construction,  mining,  agriculture,  forestry,  recreation, 
military,  tire  and  auto  companies. 

Collaboration:  Auburn,  UC  Davis,  Goodyear 

S.  Shoop 

Cold  Regions  Research  and  Engineering  Laboratory  (CRREL) 

Hanover,  NH  03755-1290 
Phone:  (603)646-4321 
Fax:  (603)646-4640 
e-mail:  sallys@hanover-crrel.army.mil 


The  Role  of  High  Resolution  Simulations  in  Vehicle  Performance  Assessment 

Roger  A  Wehage1 


Introduction 

This  technical  note  summarizes  a  proposal  previously  submitted  to  the  Marine  Corps  for 
developing  an  interactive,  high  resolution  vehicle  modeling  and  simulation  methodology  to 
support  the  acquisition  and  evaluation  of  a  future  Medium  Tactical  Vehicle  Replacement 
(MTVR)  system.  The  discussion  gives  an  idea  of  the  many  critical  subsystems  which  contribute 
to  a  vehicle’s  mobility,  ride  quality  and  dynamic  stability,  and  the  level  of  effort  required  to 
implement  and  validate  such  a  system. 

Motivation  for  High  Resolution  Interactive  Models 

A  Real-Time  Operator-ln-The-Loop  Concept 

The  System  Simulation  &  Technology  Division  (SSTD)  at  the  Tank-Automotive 
Research  &  Development  Engineering  Center  (TARDEC)  and  the  Nevada  Automotive  Test 
Center  (NATC)  proposed  to  cooperate  in  the  development  of  an  advanced  high  resolution  vehicle 
dynamic  performance  evaluation  methodology.  The  methodology  would  develop  and  implement 
procedures  to  generate  vehicle  models  capable  of  running  at  or  near  real  time  on  moderate  to 
large  scale  computers.  The  real-time  simulation  models  would  provide  rapid  and  accurate 
dynamic  performance  evaluations  because  the  critical  human  decision  making  and  response 
processes,  which  are  impossible  to  characterize  and  model,  would  be  input  directly  to  the  models 
by  a  user  through  a  graphics-based  vehicle  simulation  workstation.  The  optimized  vehicle 
equations  and  computer  algorithms  would  be  generated  by  a  new  modeling  and  simulation 
methodology  under  development  at  TARDEC,  called  Symbolically  Optimized  Vehicle  Analysis 
System  (SOVAS).  SOV AS-generated  equations  would  be  executed  under  control  of  an 
interactive  graphics-based  program  called  Dynamic  Response — Interactive  Vehicle  Emulator 
(DRIVE).  DRIVE  would  provide  a  continuous  display  of  the  driver’s  view  of  the  surrounding 
vehicle  and  terrain  operating  scenarios  and  would  allow  him  to  input  control  commands  from  a 
workstation  interface,  similar  to  those  in  the  vehicle.  All  vehicle  subsystems,  and  vehicle/terrain 
interactions  would  be  modeled  as  accurately  as  feasible  to  achieve  the  maximum  possible 
resolutions  while  still  retaining  real-time  simulation  capability.  SOVAS  and  DRIVE-based 
vehicle  models  would  provide  performance  discrimination  capabilities  which,  until  now,  could 
only  be  obtained  through  expensive,  time  consuming,  and  possibly  dangerous  field  and 
laboratory  tests.  A  second  main  feature  was  that  the  methodology  would  also  allow  application 
of  the  same  high  resolution  performance  assessments  to  concepts  and  prototypes  that  may  only 
exist  on  paper  or  that  may  be  unavailable  for  extensive  field  and  laboratory  testing. 

Many  Subsystems  are  Critical  to  Vehicle  Performance  Predictions 

A  high  performance  tactical  wheeled  vehicle  contains  many  unique  subsystems,  of  which 
a  good  number  are  critical  to  its  successful  performance.  Many  of  these  subsystems  significantly 
influence  a  vehicle’s  dynamic  performance,  and  an  operator’s  ability  to  control  it  Various 
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subsystems  are  more  critical  than  others  for  a  given  performance  aspect,  and  it  is  generally 
impossible  to  specify  their  relative  degree  of  importance.  Thus,  the  best  that  can  be  done  in  lieu 
of  a  time  consuming  and  expensive  sensitivity  study,  is  to  model  each  subsystem  to  the  highest 
degree  of  resolution  commensurate  with  the  modeling  and  simulation  capabilities.  Validation 
studies  which  encompass  the  important  dynamic  performance  aspects  are  necessary  to  insure  that 
each  subsystem  model  is  correct  and  has  the  necessary  resolution. 

Successful  vehicle  performance  is  becoming  more  and  more  dependent  on  the  application 
of  on-board  sensors  and  computer  algorithms  to  carry  out  many  of  the  critical  tasks.  Successful 
performance  prediction  models  must  also  include  accurate  representation  of  the  sensor  and 
controller  dynamics,  as  well  as  the  onboard  computer  algorithms.  Thus  the  next  generation 
modeling  and  simulation  methodologies  must  place  heavy  emphasis  on  accurate  characterization 
of  these  subsystems,  and  modeling  and  simulating  them. 

A  vehicle’s  chassis  and  the  driver/passenger  compartment  are  the  framework  of  its 
operating  system.  The  vehicle  interacts  with  its  surrounding  environment  and  the  operators 
through  a  number  of  intricate  subsystems.  Likewise  the  operators  interact  with  the  vehicle  and 
the  environment  by  sensing  and  controlling  many  of  these  subsystems.  While  it  is  feasible  to 
accurately  model  the  subsystems  in  a  graphics-based  workstation,  it  is  impossible  to  model  the 
operators.  The  most  cost  effective  alternative  is  to  allow  the  operator  to  control  a  vehicle  model 
directly  from  a  workstation.  Graphical  feedback  would  allow  the  driver  to  “see”  and 
“experience”  what  is  happening,  both  inside  and  outside  the  vehicle.  Various  controls  such  as 
joysticks,  switches  and  on-screen  touch  sensitive  devices  would  allow  the  driver  to  control  the 
vehicle.  Since  a  vehicle  model’s  ‘driver’  would  not  be  in  the  physical  environment,  he  would 
not  get  many  of  the  normal  visual,  audio  and  touch  sensory  feedback  signals.  Thus  the  DRIVE 
methodology  would  have  to  be  designed  to  augment  the  graphical  system  with  additional  visual 
and  audio  signals.  Experience  with  ground  and  flight  simulators  has  shown  that  operators  can 
effectively  control  systems  with  a  reduced  set  of  sensory  inputs. 

Critical  Elements  of  a  Vehicle  Model 

The  following  discussion  gives  a  brief  description  of  the  most  important  MTVR 
subsystems.  This  discussion  gives  a  general  indication  of  the  critical  interaction  dynamics  which 
goes  on  between  these  systems,  the  operator,  the  vehicle,  and  its  surrounding  environment  The 
subsystems  have  not  been  arranged  in  any  particular  order  of  importance.  In  fact  it  is  impossible 
to  rank  them  because  it  is  only  their  overall  synergistic  performance  that  can  be  evaluated.  It  is 
also  impossible  to  test  each  subsystem  separately  because  they  are  so  tightly  coupled  together 
and  interdependent. 

Role  of  the  Drive  Train 

The  engine  model  is  discussed  first  because  it  provides  the  driving  energy  source  for  the 
vehicle.  It  is  partially  controlled  by  the  driver  through  a  number  of  inputs,  and  by  other 
subsystems  that  may  be  partially  under  computer  control.  The  primary  control  is  through  the 
throttle  which  may  be  indirectly  affected  by  a  governor  and  other  electronic  control  devices. 
Computer  algorithms  may  directly  or  indirectly  control  a  number  of  engine  parameters  based  on 
various  sensed  vehicle  states.  The  driver  may  also  select  various  levels  of  engine  braking 
involving  jake  brakes  or  other  types  of  installed  engine  braking  devices.  Engine  braking  may 
also  be  influenced  by  the  state  of  operator  applied  brakes  and  the  antilock  braking  system  (ABS). 
It  is  not  feasible  to  develop  and  incorporate  high  resolution  dynamic  engine  models  into  a  real¬ 
time  simulation  methodology,  so  engine  torque  output  should  be  based  on  empirically  measured 
data  curves  which  depend  on  a  number  of  controller  and  vehicle  state  variables. 
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Most  of  the  systems  in  a  vehicle  rely  on  hydraulic,  air  and  electrical  supply  systems  to 
function  properly.  These  systems  must  be  characterized  and  modeled  to  insure  they  are  capable 
of  providing  the  necessary  energy  to  drive  all  of  the  actuators.  The  power  supply  models  will 
determine  the  power  drain  on  the  engine. 

A  rotating  inertia  would  be  used  to  represent  the  engine  crankshaft  and  flywheel  which  is 
connected  to  the  transmission  torque  converter.  A  typical  torque  converter  assembly  consists  of 
a  rotary  impeller,  turbine  and  stator  designed  to  smooth  out  transient  motions  in  the  drive  train 
and  provide  the  torque  gains  necessary  to  rapidly  bring  a  vehicle  up  to  speed  when  accelerating, 
or  to  slow  it  down  during  engine  braking.  An  electro-hydraulically  controlled  clutch  is  used  to 
lock  out  the  torque  converter  when  it  saturates,  to  improve  engine  efficiency.  As  noted  above  for 
the  engine,  it  is  not  feasible  to  develop  and  incorporate  high  resolution  dynamic  torque  converter 
models  into  a  real-time  simulation  methodology,  so  the  turbine  and  clutch  output  ton]ue 
relationships  are  best  determined  empirically  and  incorporated  into  the  model  as  functions  of 
various  controller  and  system  states. 

The  torque  converter  output  is  coupled  into  the  vehicle’s  transmission  gear  box.  A 
number  of  user-selectable  or  computer  controlled  gear  reduction  systems,  which  are  activated  or 
deactivated  by  computer  controlled  electro-hydraulic  clutches,  work  to  give  the  vehicle  a  suitable 
speed  and  power  range.  Again,  it  is  not  feasible  to  develop  and  incorporate  high  resolution 
clutch  torque  models  into  a  real-time  simulation  methodology,  so  these  outputs  are  best 
empirically  determined  and  included  in  the  model  along  with  models  of  the  control  computer 
algorithms  (or  the  actual  control  computer  code  from  the  vehicle)  and  other  measured  system 
states. 


The  MTVR  has  an  intricate  arrangement  of  inter  and  intra  axle  differentials  to  give  a 
wide  range  of  possible  axle  drive  configurations.  The  differentials  are  electro-pneumatically 
activated  and  deactivated  through  manual  switches  located  in  the  cab.  Time  constants  for  the 
actuators  would  be  determined  empirically.  The  driver  would  be  able  to  control  the  state  of  these 
systems  according  to  his  observance  of  the  surrounding  environment  and  the  vehicle’s  dynamic 
states. 


The  MTVR  has  a  unique  suspension  system,  as  do  other  high  mobility  vehicles.  They 
represent  intricate  interactions  between  axles,  springing  devices,  dampers,  wheels,  travel  limiters, 
tires  and  the  terrain.  In  fact,  the  suspension  must  be  considered  an  integral  part  of  the  power 
train  from  the  engine  all  the  way  to  the  tires  and  ground.  Poorly  designed  suspension  systems 
(and  models)  could  cause  internal  oscillations  resulting  in  loss  of  traction  control,  poor  ride 
quality,  excessive  wear,  etc.  Thus  it  is  important  to  accurately  model  the  suspensions,  as  well  as 
the  other  interconnected  systems,  to  insure  that  the  model,  itself,  does  not  inadvertently  degrade 
or  improve  a  vehicle’s  predicted  performance. 

The  steering  system  also  affects  a  vehicle’s  dynamic  performance.  It  is  necessary  to 
accurately  represent  the  steering  kinematics  all  the  way  from  the  steering  wheel  to  the  wheel 
hubs.  Incorrect  steering  kinematics  and  compliance  can  result  in  over  steering,  under  steering  or 
incorrect  load  transfer  to  the  suspensions  and  tires.  This  can  cause  vehicle  instabilities,  loss  of 
traction  and  other  problems.  The  time  response  haracteristics  of  hydraulically  power  assisted 
units  will  also  affect  steering  performance. 

Between  the  Vehicle  and  the  Road 

The  MTVR  central  tire  inflation  system  (CTIS)  has  been  designed  to  provide  a  means  for 
automatically  adjusting  tire  pressures  in  response  to  current  terrain  and  operating  conditions. 
Depending  on  the  system  design,  the  CTIS  states  may  be  manually  controlled  from  the  cab,  or 
they  may  be  directly  or  indirectly  controlled  by  one  of  the  on-board  computers  used  for  other 
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purposes.  The  response  characteristics  of  the  CTIS  control  valves  would  be  determined 
empirically. 

On  the  MTVR,  ABS  is  designed  to  control  wheel  slip  during  manual  braking.  It  senses 
wheel  rotational  speeds  and  brake  pad  temperatures,  and  uses  this  information  to  control  the 
pneumatically  actuated  brakes.  Some  systems  may  also  use  ABS-type  control  algorithms  to 
improve  traction  control  while  accelerating  or  towing  loads.  ABS  requires  a  sophisticated 
computer  program  to  perform  its  functions,  and  this  program  would  have  to  be  included  in  the 
vehicle  model  as  well. 

A  tire  is  one  of  the  most  complex  and  critical  subsystems  in  a  vehicle  and  a 
comprehensive  representation  of  its  interaction  dynamics  with  irregular  nondeforming  or 
deformim  aces  has  not  been  obtained,  even  through  extensive  research  efforts.  Thus  the 
most  feai  solution  is  to  use  models  based  on  empirically  determined  data.  A  dynamic 
enveloping  ure  model  could  be  used  to  obtain  a  set  of  tire  state  variables  which  would  then  be 
inserted  into  an  empirical  tire/soil  interaction  model.  This  model  would  return  a  set  of 
interaction  forces  and  moments  which  would  then  be  applied  back  to  the  dynamic  enveloping  tire 
model.  An  iterative  procedure  would  be  used  to  converge  the  two  models  to  a  common  solution. 
This  approach  would  require  very  low  order  models  to  obtain  accurate  results  so  real-time 
simulations  would  still  be  possible. 

Accurate  representation  of  tire/soil  interaction  dynamics  is  very  crucial  to  many  of  the 
dynamic  vehicle  performance  analyses.  Mobility  is  likely  the  most  important  vehicle 
performance  aspect  and  many  of  the  major  vehicle  subsystems  were  designed  with  this  in  mind. 
Dynamic  stability  and  ride  quality  are  two  other  important  performance  aspects  which  depend  on 
many  of  the  subsystems,  and  tiie/soil  interaction  dynamics.  Thus  it  is  clear  that  if  computer- 
based  analysis  is  to  be  a  useful  tool  for  studying  these  phenomena,  the  corresponding  models 
must  be  accurate  enough  to  emulate  the  critical  subsystems’  dynamic  response  characteristics. 

Empirical  characterization  of  tire/soil/terrain  interaction  dynamics  is  encumbered  by  the 
wide  range  of  tire,  soil  and  terrain  configurations.  A  program  to  accomplish  such  a  task  would 
require  identification  of  the  most  important  tire,  soil  and  terrain  parameters  and  elimination  of  all 
others.  Then  a  comprehensive  parameter  identification  and  measurement  procedure  would  have 
to  be  established. 

Interfacing  a  Vehicle  Model,  the  Driver  and  the  Environment 
The  Role  of  a  Vehicle  Operator  in  Performance  Prediction 

This  discussion  has  emphasized  the  importance  of  accurate  representation  of  the  various 
vehicle  subsystem  models.  However,  the  most  important  vehicle  subsystem  which  is  impossible 
to  characterize  or  model  is  the  driver.  The  best  driver  model  which  might  possibly  be  defined 
could  severely  compromise  a  vehicle’s  performance  because  it  simply  could  not  come  close  to 
emulating  a  human’s  thought  and  response  processes.  The  most  feasible  solution,  which  would 
allow  a  vehicle  model  to  be  accurately  exercised  through  the  widest  range  of  highly  nonlinear 
displacements  that  a  corresponding  vehicle  in  the  field  might  experience,  is  to  allow  an  operator 
to  perform  these  functions,  in  real  time,  at  a  graphics  workstation  or  in  some  other  physical 
simulation  environment. 

Interfacing  the  Operator  With  the  Vehicle 

Two  major  efforts  would  be  required  for  interfacing  the  driver  to  the  vehicle  simulation 
model.  First,  adequate  controls  would  have  to  be  provided  to  allow  him  to  comfortably  and 
accurately  input  the  necessary  commands  to  the  vehicle.  Second,  sufficient  graphical  and  audio 
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feedback  signals  would  have  to  be  provided  so  he  could  continuously  determine  what  the  vehicle 
is  doing.  Numerous  ground  vehicle  and  flight  simulators  have  demonstrated  that  it  is  possible  to 
provide  sufficient  controls  and  sensory  inputs,  and  that  a  user  can  adequately  control  such 
systems  in  an  emulated  environment  A  number  of  field  and  laboratory  tests  with  actual  vehicles 
would  have  to  be  performed  to  leam  what  controls  and  sensory  inputs  are  necessary  for  adequate 
vehicle  control. 

Interfacing  the  Operator  With  the  Environment 

In  a  realistic  operational  environment,  it  would  also  be  necessary  to  account  for  other 
types  of  vehicle  and  operator  interactions  with  the  surrounding  environment.  The  most 
important  interaction  would  be  chassis  interference  with  the  terrain  and  surrounding  vegetation. 
An  extensive  amount  of  empirical  data  has  been  collected  for  the  NATO  Reference  Mobility 
Model  (NRMM)  which  could  be  used  to  support  such  a  modeling  and  simulation  effort  In 
addition,  some  commercial  software  programs  could  be  used  to  incorporate  Defense  Mapping 
Agency  Digital  Terrain  Elevation  Data  (DMA-DTED)  and  Digital  Feature  Analysis  Data  (DMA- 
DFAD)  directly  into  MTVR  models.  NAT C  has  access  to  world-wide  DTED  and  DFAD  data 
bases  and  procedures  could  be  developed  to  extend  their  resolutions  down  to  any  degree  of 
accuracy  necessary  to  carry  out  the  corresponding  analyses.  DTED  data  files  give  the  terrain 
elevation  data  for  many  locations  in  the  world  and  any  number  of  surface  and  material  properties 
could  be  associated  with  each  piece  of  the  terrain.  DFAD  data  files  define  the  types  of  features 
on  the  terrain  profiles  and  their  properties.  For  example,  DFAD  files  have  attributes  defining 
feature  type,  height,  orientation,  identification,  surface  material  code,  etc.,  and  any  number  of 
additional  feature  properties  could  be  added  as  necessary  to  extend  this  data  base.  A  number  of 
dynamic  form  drag  and  interference  models  have  been  reported  in  the  literature.  These  models 
could  be  investigated  and  incorporated  into  the  DRIVE  methodology  to  obtain  representative 
vehicle/obstacle/vegetation  interference  models.  This  data  would  be  important  for  providing  the 
necessary  vision  impairments  to  vehicle  operators  as  well. 

Importance  of  Trailer  and  Wagon  Models 

Equal  emphasis  must  be  placed  on  the  characterization  of  trailers  and  wagons,  and  the 
development  of  high  resolution  models  for  these  systems.  Trailer  dynamics  can  adversely  affect 
a  truck’s  mobility,  ride  quality  and  dynamic  stability.  In  fact,  any  truck/trailer  combination  must 
be  considered  as  a  single  integrated  system  with  its  own  mobility,  ride  quality  and  dynamic 
stability  properties.  Failure  of  the  trailer  generally  means  failure  of  the  system  as  well.  If  a 
given  trailer  uses  the  truck’s  hydraulic,  air  or  electrical  supply  systems,  then  these  extra  loads 
must  also  be  accounted  for  in  the  truck  models.  If  any  trailer  functions  are  influenced  or 
controlled  by  the  truck’s  on  board  computers,  such  as  ABS  functions,  procedures  must  be 
established  to  incorporate  them  into  the  models  as  well.  Finally,  procedures  must  be  established 
to  give  the  driver  adequate  visual  and  sensory  feedback  on  the  trailer’s  dynamic  state  so  he  can 
respond  accordingly. 

System  Characterization  and  Model  Validation 

System  characterization  and  model  validation  is  considered  an  extremely  important  and 
integral  part  of  every  subsystem  model  development  effort.  TARDEC  engineers  have  an 
extensive  background  in  defining  the  supporting  mathematical  equations,  developing  the 
computer  algorithms,  obtaining  the  system  states  and  interpreting  the  results.  However,  they 
have  very  limited,  or  no  knowledge  of  the  functional  operation  of  most  of  the  MTVR’s  complex 
subsystems.  On  the  other  hand,  NATC  engineers  have  an  extensive  knowledge  of  the  intricate 
operation  of  every  major  subsystem  because  they  were  either  responsible  for  the  system 
performance  specifications  to  other  developers,  or  they  developed  the  systems  themselves.  They 
also  have  close  working  relationships  with  many  of  the  vendors  and  system  developers.  For  a 
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proposed  program  to  be  successful,  TARDEC  engineers  would  have  to  leam  the  functional 
operation  of  these  subsystems  and  NATC  engineers  would  have  to  leam  the  basics  of  computer- 
based  vehicle  modeling  and  simulation  so  the  two  groups  could  create  the  best  possible  vehicle 
models  and  analysis  tools.  Thus  TARDEC  and  NATC  engineers  would  have  to  work  closely  at 
each  step  of  a  cooperative  project  to  insure  that  the  operational  performance  of  every  major 
subsystem  is  thoroughly  understood  by  everyone  involved.  In  addition,  the  participants  would 
have  to  insure  that  all  necessary  subsystem  characterization  and  model  validation  tests  were 
carefully  defined  and  performed  so  the  models  could  be  accurately  defined  and  validated. 
TARDEC  and  NATC  engineers  would  have  to  leam  from  each  other  in  such  a  project  so  they 
would  be  in  a  good  position  to  apply  this  gained  knowledge  and  experience  to  follow-on  efforts. 

Motivation  for  Support  of  the  Acquisition  Process 

It  was  expected  that  such  an  effort  would  result  in  a  comprehensive  performance 
specification  which  would  go  into  a  subsequent  Request  For  Proposal  (RFP)  and  that  the  RFP 
would  require  each  prospective  bidder  to  extensively  evaluate  his  concepts  using  high  resolution 
computer-based  vehicle  models.  It  was  further  expected  that  this  effort  would  give  NATC 
engineers  the  ability  to  define  and  execute  these  high  resolution  models  (with  some  assistance 
from  TARDEC  engineers),  and  that  most  serious  contractors  would  seek  help  from  NATC.  It 
was  also  expected  that  TARDEC  engineers  would  be  tasked  to  model,  simulate  and  evaluate  the 
proposed  systems  following  receipt  of  the  proposals,  and  during  further  down-select  phases.  The 
critical  issue  for  TARDEC’ s  evaluation  efforts  has  always  been  the  lead  time  and  manpower 
required  to  develop  and  validate  comprehensive  models  in  order  to  give  reasonable  proposal 
evaluation  response  times.  If  contractors  were  encouraged  to  use  such  a  methodology  in  the 
initial  development  phases,  and  if  NATC  engineers  were  to  use  their  gained  knowledge  to  insure 
that  the  contractors  carefully  obtained  the  necessary  data  for  the  models,  then  they  would  be  in 
an  excellent  position  to  supply  accurate  model  data  or  entire  models  to  TARDEC  in  a  timely 
manner.  This  would  have  the  added  advantage  of  forcing  the  contractors  to  create  better 
prototype  designs.  In  addition,  it  would  allow  TARDEC  engineers  to  evaluate  the  proposals  in  a 
more  timely  manner,  thereby  increasing  the  probability  the  Marine  Corps  would  get  the  best 
possible  product. 

Development  and  Validation  of  the  Methodology — A  Cooperative  Effort 

The  second  major  task  would  concentrate  on  the  development  and  validation  of  a 
comprehensive  high  resolution  model  of  the  MTVR  and  trailers.  In  this  process,  the 
methodologies  to  support  high  resolution,  real-time,  man-in-the-loop  simulations  in  a  graphics- 
based  workstation  environment  would  be  perfected.  The  first  effort  would  concentrate  on 
expanding  the  modeling  and  simulation  methodology  in  preparation  for  defining  the  performance 
specifications  and  evaluation  criteria  which  would  go  into  the  Request  For  Proposal.  Six  major 
research  areas  were  anticipated  to  complete  this  effort: 

1 .  Validate  derivability  of  the  methodology  against  extensive  MTVR  field  test  data. 

2.  Expand  the  methodology  to  assess  MTVR  performance  in  the  Marine  Corps 
world-wide  operational  scenarios. 

3.  Define  and  develop  a  unique  set  of  computer-based  vehicle  simulations  which 
must  be  carried  out  to  quantify  and  discriminate  between  new  concept  and  prototype  vehicle 
performance  parameters. 

4.  Develop  advanced  subsystem  characterization,  model  development  and  validation 
procedures  necessary  to  minimize  the  time  between  receipt  of  vehicle  data  and  completion  of 
performance  analyses. 
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5.  Develop  comprehensive  vehicle  performance  specifications,  data  requirements 
and  evaluation  procedures  Data  Item  Descriptions  (DIDS)  which  must  be  incorporated  into  the 
RFP’s  to  insure  that  sufficient  and  accurate  information  is  provided  to  and  by  contractors. 

6.  Design  flexibility  into  the  modeling  and  simulation  methodology  to  simplify  the 
definition  and  incorporation  of  all  unique  design  features  of  each  proposed  concept  or  prototype 
into  the  models,  thus  minimizing  the  possibility  of  compromising  promising  new  design  features. 

Field  Tests  for  Model  Validation 

The  above  discussion  clearly  indicates  that  the  SOVAS-generated  MTVR  and  trailer 
models  would  have  to  be  quite  extensive  and  sophisticated  in  order  to  address  all  of  the  critical 
dynamic  operational  performance  aspects  of  a  system  of  this  magnitude.  Furthermore,  the 
DRIVE-based  man-in-the-loop  methodology  would  have  to  be  capable  of  insuring  that  the 
models  can  be  accurately  exercised  to  answer  the  critical  design  and  evaluation  questions.  Thus 
the  first  and  most  critical  step  of  this  effort  would  be  to  defme  and  carry  out  a  battery  of 
comprehensive  and  carefully  controlled  field  tests  to  generate  SOVAS  and  DRIVE  validation 
data.  The  scope  of  these  tests  would  be  determined  in  conjunction  with  the  actual  model 
development  efforts.  In  general,  it  would  be  impossible  to  predict  the  type  and  detail  of  model 
validation  tests  required  until  the  subsystem  characterization  and  model  development  efforts 
wert  underway  or  completed.  Thus  part  of  the  effort  would  be  to  identify  and  define  the 
validation  procedures  that  would  have  to  be  conducted.  The  confidence  gained  in  the  subsystem 
characterization  and  model  development  procedures  from  these  validation  efforts  would  allow 
TARDEC  and  NATC  engineers  to  reliably  extend  the  methodologies  to  new  concepts  and 
prototypes,  and  to  new  operational  scenarios  as  described  above. 

The  Marine  Corps  world-wide  MTVR  operational  performance  requirements  involve  a 
large  number  of  terrain,  soil,  vegetation  and  obstacle  types.  NATC  and  TARDEC  would  create 
extensive  DTED  and  DFAD  data  bases  for  these  regions.  MTVR  operational  performance  on 
these  data  bases  would  be  validated  against  test  data  within  these  regions  or  test  data  taken  from 
other  areas  with  similar  types  of  terrain,  soil,  vegetation  and  obstacles.  The  validated  model  and 
data  bases  would  then  be  available  for  future  concept  and  prototype  evaluations. 

NATC  and  other  Government  agencies  perform  a  broad  spectrum  of  carefully  controlled 
field  tests  in  order  to  quantify  a  vehicle’s  dynamic  operational  performance  capabilities.  These 
tests  are  run  on  many  types  of  terrain  profiles,  conditions  and  operating  environments.  All  field 
tests  which  could  be  augmented  or  replaced  by  SOVAS  and  DRIVE-based  vehicle  emulations 
should  be  defined  and  validated.  Test  procedures  would  have  to  be  identified  and  installed  into 
the  DRIVE  workstation  environment  It  would  be  best  if  all  such  procedures  were  in  place  prior 
to  receipt  of  proposals  and  prototypes  in  order  to  achieve  the  necessary  evaluation  response 
times. 

System  Characterization,  Model  Development  and  Validation 

Another  critical  area  that  would  have  to  be  addressed  by  TARDEC  and  NATC  engineers 
is  the  definition  and  development  of  advanced  subsystem  characterization,  model  development 
and  validation  procedures  to  minimize  the  time  between  receipt  of  vehicle  data  and  completion 
of  performance  analyses.  It  would  be  ideal  if  every  contractor  provided  data  or  subsystem 
models  which  could  be  entered  directly  into  the  SOVAS  and  DRIVE  methodologies.  However, 
TARDEC’ s  experience  is  that  many  contractors  either  lack  the  ability  or  desire  to  provide  data  or 
models,  and  when  such  information  and  models  are  supplied,  their  accuracy  and  resolution  is 
often  questionable.  Also,  contractors  often  consider  their  on-board  control  algorithms  and 
systems  to  be  proprietary  and  would  refuse  to  release  information  to  TARDEC  or  NATC.  If  high 
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resolution  vehicle  simulations  are  to  be  useful  evaluation  tools,  the  RFP’s  must  be  very  specific 
in  what  the  contractors  must  supply  and  how  the  information  would  be  used  in  the  evaluation 
process.  TARDEC  and  NATC  engineers  would  have  to  clearly  define  what  is  required  to 
support  these  modeling  and  simulation  efforts  and  what  the  contractors’  responsibilities  would 
be. 

Development  for  Support  of  the  Acquisition  Process 
Defining  the  Performance  Specifications 

TARDEC  and  NATC  engineers  would  also  have  to  develop  comprehensive  vehicle 
performance  specifications,  data  requirements  and  evaluation  procedures  DEDS  which  must  be 
incorporated  into  the  RFP’s  to  insure  that  sufficient  and  accurate  information  is  provided  to  the 
contractors.  An  important  goal  of  this  proposed  program  would  be  to  quantify  unique  vehicle 
design  features  that  contribute  to  improved  mobility,  stability  and  ride  quality,  and  to  develop 
computer-based  modeling  and  simulation  methodologies  that  would  allow  rapid  evaluation  of 
new  concepts  and  prototypes  to  determine  how  well  the  contractors’  proposals  have  met  these 
criteria.  The  challenge  would  be  the  wording  of  these  specifications  so  the  Government  would 
not  be  telling  the  contractors  how  to  design  their  vehicles,  yet  would  have  a  good  probability  that 
the  proposals  would  meet  the  Marine  Corps  needs.  The  DIDS  would  have  to  be  very  specific  on 
what  data  the  contractors  must  supply  and  in  what  forms  it  must  be  supplied.  Furthermore,  it 
should  be  made  very  clear  that  incomplete  or  inaccurate  data  could  lead  to  incorrect  evaluation  of 
their  proposals,  and  possible  rejection. 

Model  Adaptability  to  Unique  Vehicle  Designs 

The  proposed  effort  would  also  design  the  SOVAS  and  DRIVE  methodologies  to  accept 
a  wide  range  of  subsystem  models.  This  would  be  necessary  in  order  to  simplify  the  definition 
and  incorporation  of  all  unique  design  features  of  each  proposed  concept  or  prototype  into  the 
models,  thus  minimizing  the  possibility  of  compromising  promising  new  design  features.  No 
features  of  the  MTVR  or  trailer  models  would  be  hard  coded  into  the  SOVAS  and  DRIVE 
methodologies.  Each  of  the  major  subsystems  which  the  MTVR,  and  nearly  every  concept  or 
prototype  must  have,  would  be  defined  as  stand-alone  modules  which  could  be  tailored  and 
optimized  for  each  vehicle  model.  That  is,  each  major  subsystem  would  have  a  generic  module 
with  the  main  functional  requirements  which  could  be  tailored  and  optimized  as  necessary.  A 
major  goal  of  this  project  would  be  to  gain  experience  with  the  subsystem  characterization, 
model  development  and  validation  efforts  using  the  MTVR  so  it  would  become  more  of  a  routine 
procedure  to  tailor  and  optimize  the  generic  modules  in  the  follow-on  evaluation  efforts. 
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Soil  Plowing  Using  die  Discrete  Element  Method 


David  A.  Homer 

U.S.  Army  Engineer  Waterways  Experiment  Station 
Vicksburg,  Mississippi  39180-6199 


BACKGROUND 

Many  vehicle-soil  interaction  problems  involve  large  discontinuous 
deformations.  Such  problems  include  sinkage  of  vehicle  tracks  and  tires  in  soft 
soils  and  vehicle  plowing  problems.  The  traditional  use  of  continuum  mechanics  is 
limited  to  large  strains  with  continuous  deformation.  The  principal  difficulty 
comes  from  the  mathematical  description  of  the  kinematics  that  describes  the 
movement  of  material  "particles"  within  die  continuum.  In  a  continuum, 
movement  must  obey  compatibility  relationships  that  preclude  formation  of  slip 
planes  and  separations.  However,  in  real  materials,  large  deformations  can  occur 
that  violate  compatibility  in  the  strict  sense  imposed  by  continuum  formulations. 

To  model  discontinuous  deformations,  as  done  in  finite  element  analysis  of  fracture 
propagation  or  at  interfaces  in  the  soil-structure  interaction  problems,  special  joint 
elements  are  often  introduced  to  model  slip  or  separation  (Wong  and  Hanna,  1977). 

Use  of  these  elements  complicate  analyses  for  small  deformation  problems  and 
become  excessively  complicated  for  large  deformation  problems.  The  joint 
elements  also  require  assumptions  to  be  made  on  the  location  of  the  discontinuities. 

An  alternative  to  the  continuum  description  for  soil  and  rock  mechanics 
problems  is  the  Discrete  Element  Method  (DEM)  which  models  the  material  as  a 
collection  of  individual  unconnected  particles.  The  motion  of  the  particles  is 
controlled  by  Newton's  laws  of  motion  and  is  only  restricted  kinematically  by  the 
requirement  that  particles  cannot  penetrate  each  other.  Several  authors  have  used 
DEM  to  model  granular  assemblies  (Cundall  and  Strack,  1979,  Christoffersen,  et. 
Al.,  1981,  Bathurst  and  Rothenburg,  1988).  Ng  and  Dobiy  (1991)  used  DEM  to 
model  small  strain  cyclic  loading.  Their  simulation  results  agreed  closely  with 
trends  found  in  laboratory  tests  of  sands.  Shukla  and  Sadd  (1990)  used  DEM  to 
investigate  how  mechanical  stress  waves  propagate  in  granular  material  and  how 
they  are  influenced  by  media  microstructure.  DEM  has  been  used  to  model  the 
results  of  tunnel  failure  resulting  of  a  nuclear  explosion  (Heuze,  et.  Al.,  1991). 
Sophisticated  algorithms  have  been  developed  to  describe  the  evolution  of  the 


particulate  system  including  the  formation  and  breaking  of  inter-element  contacts 
(Ting,  et.  Al.,  1989).  Computing  forces  between  elements  requires  relationships  to 
describe  normal  and  shear  interaction  at  the  contacts.  In  some  models  the 
individual  particles  can  "break"  when  stress  conditions  within  the  particle  reach 
some  critical  level  (Cundall  and  Hart,  1985).  Typically  soil  particles  have  been 
modeled  as  two-dimensional  circular  rigid  disks.  Ting  (1991)  has  developed  an 
ellipse-based  two-dimensional  particle  to  represent  contact  flatness  and  particle 
angularity.  Six-sided  solid  shapes  have  been  used  to  model  granular  material 
(Ghaboussi  and  Barbosa,  1990).  The  predominant  disadvantage  of  DEM  is  the 
enormous  computational  requirement,  which  is  a  result  of  keeping  track  of  all 
particle  contact  locations. 

PARTICLE  PLOW  SIMULATION 

A  1,200  particle  simulation  was  performed  to  illustrate  the  potential  use  of 
the  DEM  for  investigating  the  response  of  the  particles  to  a  rigid  plow  moving 
through  the  particulate  mass  .  The  simulation  consisted  of  dropping  the  particle 
into  the  test  chamber,  allowing  the  particle  come  to  rest,  and  then  moving  the  plow 
through  the  test  specimen  at  a  constant  rate.  The  plow  blade  angle  was  set  to  45°. 
Figure  1  plots  the  force  magnitude  of  the  particles  prior  to  plowing.  The  lighter 
the  color  the  more  force  acting  on  the  particle.  Figure  2  plots  the  force  magnitude 
of  the  particles  after  plowing.  Particles  near  the  plow  tip  show  an  build  up  of 
force  concentration. 

Figure  3  shows  only  those  particles  that  have  horizontal  motion  at  the  end 
of  the  plowing  test.  The  lighter  color  represents  those  particles  with  horizontal 
velocity.  Using  scientific  visualization  to  animate  the  results  of  the  simulation 
allows  for  the  observation  of  soil  deformation  patterns  forming  and  identify 
formation  of  failure  surfaces. 


CONCLUSIONS 

The  Distinct  Element  Method  is  an  analytical  tool  for  fundamental  research 
into  the  behavior  of  granular  material.  The  DEM  allows  the  simulation  of 
complex  nonlinear  interaction  problems  in  terramechanics.  The  major  drawback  to 
the  DEM  is  the  enormous  computional  requirement.  Future  research  in  the  use  of 
scaling  principle  for  full  scale  modeling  and  the  use  massively  scalable  parallel 
computer  systems  will  be  attempts  to  address  this  problem.  The  advantage  of  the 
DEM  is  that  slip  planes  and  separations  can  form  between  groups  of  particles  thus 
capturing  evolving  failure  mechanisms  in  a  simpler  and  more  realistic  way  than 
models  based  on  a  continuum  description  of  the  soil  mass. 


Figure  3.  Particle  Forces  After  Plowing 
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